Root canal feature form extraction method, device and equipment based on CBCT data and medium

Abstract

The embodiment of the invention discloses a root canal feature form extraction method, device and equipment based on CBCT data and a medium, and belongs to the technical field of stomatology and medical image.The method comprises the steps that filtering denoising and dental pulp segmentation processing are carried out on a CBCT image, and surface drawing three-dimensional reconstruction is carried out on a tooth root canal to generate a three-dimensional model; based on regular point cloud data in the three-dimensional model, adopting a point-surface method to realize curved surface reconstruction, and generating a root canal entity model; and performing coordinate alignment processing on the root canal solid model in a three-dimensional space, and measuring the length of the root canal, the direction of a root canal orifice, the curvature of the root canal and the depth of dental pulp by fitting the approximate center line of the root canal. According to the method, the corresponding three-dimensional model is quickly reconstructed based on the CBCT cross-sectional image obtained by scanning the teeth of the patient, and related data such as the dental pulp depth and the size of the root canal orifice are acquired from the three-dimensional model of the teeth so as to assist the oral surgery.

Description

Technical Field

[0001] The present application relates to a method, device, equipment and medium for extracting characteristic morphology of root canals based on CBCT data, and belongs to the field of stomatology and medical imaging technology. Background Art

[0002] Accurately measuring the morphological data of the dental pulp is an important part of clinical oral surgery. The tooth structure varies with age, gender, race, etc., especially the structure of the dental pulp and root canal is very complex. These tooth morphological data can not only provide doctors with risk assessment and efficient surgical plans before surgery, such as establishing the optimal pulp opening position and direction, determining the risk of dentin fracture, and locating the root canal opening; they can also provide feedback on the surgical progress and related information during the surgery, such as the progression of the patient's condition, visual information of the surgical visualization, and mechanical analysis of the instrument working in the dental pulp.

[0003] CBCT data is used to display the entire tooth 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 projection: CBCT uses KV-level X-rays for projection. The radiation source adopts three-dimensional cone beam X-rays to perform circular digital projection around the projection object.

[0005] Step 2, data acquisition: collect two-dimensional projection data at multiple angles through a flat-panel detector.

[0006] Step 3, data reconstruction: reconstruct the collected two-dimensional projection data through a computer, and use a cone beam CT reconstruction algorithm (such as the FDK algorithm) to reconstruct the two-dimensional projection data into a three-dimensional image.

[0007] The CBCT data acquired by existing equipment cannot intuitively reflect the characteristic morphology of the root canal. Therefore, how to extract the characteristic morphology of the root canal based on CBCT data to assist oral surgery is a technical problem that needs to be solved urgently. Summary of the invention

[0008] To solve the above technical problems, the embodiments of the present application respectively provide a method, device, equipment and medium for extracting root canal feature morphology based on CBCT data, so as to obtain relevant data such as pulp depth and root canal opening size from the three-dimensional model of the tooth, which can help doctors determine the pulp opening plan.

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

[0010] According to one aspect of the embodiments of the present application, a method for extracting the characteristic morphology of dental root canals based on CBCT data is provided. The method includes:

[0011] Perform filtering and denoising processing on the CBCT image and pulp segmentation processing, and perform surface rendering three-dimensional reconstruction on the dental root canal to generate a three-dimensional model;

[0012] Based on the regular point cloud data in the three-dimensional model, use the point-plane method to achieve surface reconstruction and generate a root canal solid model;

[0013] Perform coordinate alignment processing on the root canal solid model in three-dimensional space, and use the approximate center line fitting the root canal to measure the morphological characteristics of the dental root canal; wherein, the morphological characteristics of the dental root canal include at least one of root canal length, root canal orifice direction, root canal curvature, and pulp depth.

[0014] Further, performing filtering and denoising processing on the CBCT image and pulp segmentation processing, and performing surface rendering three-dimensional reconstruction on the dental root canal to generate a three-dimensional model includes:

[0015] Based on the obtained CBCT image, use the neighborhood edge model to improve the denoising effect of median filtering, 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;

[0016] Render the surface model of the root canal from the preprocessed CBCT sequence image by surface rendering.

[0017] Further, based on the obtained CBCT image, using the neighborhood edge model to improve the denoising effect of median filtering includes:

[0018] 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:

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

[0020] In the formula, M i represents the i-th neighboring point, i is a positive integer not greater than 8, and the 8 neighboring points are respectively located above, below, left, right of the central pixel point and in the 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;

[0021] Set the threshold size to T1, and record the number of neighboring points with gray levels similar to the central pixel as n. If the pixel difference d between the central pixel and the neighboring points is less than T1, then determine that the neighboring points have gray levels similar to the central pixel; let n = n + 1. After traversing all neighboring points, if n min ≤n≤n max , b min is the minimum value of the number of neighboring points, and n max is the maximum value of the number of neighboring points, then determine that the number of neighboring 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 of all pixel points in the neighborhood. Calculate the median 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 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 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 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 > 0, then do not replace the gray value of the central pixel point. If not satisfied, replace the gray value of the central pixel point with the median value F med 。

[0029] Furthermore, based on the intercept histogram of the reciprocal cross-entropy, calculate the optimal threshold for segmenting the dental pulp, and obtain the preprocessed CBCT sequence images, including:

[0030] Based on the denoised CBCT image, use gamma transformation to enhance the image contrast, and calculate the neighborhood average image and the gradient image respectively;

[0031] Combine the neighborhood average image and the gradient image with the gray image and add them pixel by pixel to form a mixed image;

[0032] Calculate the optimal threshold for segmenting the dental pulp using the intercept histogram of the mixed image;

[0033] Compare each pixel value in the mixed image with the optimal threshold and convert it into a binary image;

[0034] Clear all unnecessary connected tissues on the image boundary of the binary image, and only retain the dental pulp area to obtain the preprocessed CBCT sequence images;

[0035] Among them, calculating the optimal threshold for segmenting the dental pulp using the intercept histogram of the mixed image includes:

[0036] Generate a three-dimensional space with gray information, neighborhood average gray information, and gradient compensation. The three-dimensional space is divided into a target region O and a background region B by a plane;

[0037] Based on the three-dimensional space, calculate the final gradient through the following formula:

[0038]

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

[0040] Represent the mixed image as:

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

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

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

[0044] Calculate the reciprocal cross entropy of the target region and the background region through the following formula:

[0045]

[0046] Wherein, E(O,B) represents the reciprocal cross entropy; P(k) represents the prior probability that the neighborhood average gray value of the mixed image is k, represents the gray mean value of the target region, represents the gray mean value of the background region;

[0047] Based on the calculated reciprocal cross entropy E(O,B), determine the optimal threshold for pulp segmentation through the following formula:

[0048]

[0049] Wherein, is the optimal threshold for pulp segmentation, and arg min is the operation of taking the index of the minimum element of the array.

[0050] Furthermore, draw the surface model of the root canal from the preprocessed CBCT sequence images by surface rendering, including:

[0051] When uniformly sampling in the x, y, and z directions of the upper and lower jaw regions in three-dimensional space with a sampling interval of Δx, Δy, and Δz, represent the volume data through a ternary function; wherein, a cube region composed of eight adjacent sampling points is a voxel;

[0052] Extract the isosurface for each voxel; wherein, the isosurface is a surface composed of points with the same attributes in three-dimensional space;

[0053] Traverse the entire volume data to find the voxels containing the isosurface, and determine the intersection positions P(x, y, z) of the isosurface with each edge of the voxel through the following formula:

[0054]

[0055] In the formula, T is the isosurface threshold, M1 and M2 are the gray values of the first vertex and the second vertex 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 intersection position coordinates of the isosurface with each edge of the voxel;

[0056] Take the triangular patches formed by the intersection points P(x, y, z) of the isosurface with each edge of the voxel, 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 the isosurface, and combine the individual isosurfaces to form an isosurface triangular network;

[0057] In order to render a better 3D effect for the isosurface triangular network, it is also necessary to select a lighting model according to the normal vectors of each triangular patch to render the isosurface triangular network to obtain the surface model of the root canal; among them, the calculation formula for the normal vectors of each triangular patch is:

[0058]

[0059] In the formula, G x , G y , G z respectively represent the gradients in the X, Y, and Z axis directions at the vertices of the isosurface,, M (x+a,y,z) , M (x-a,y,z) represents the gray values of two points (x + a, y, z) and (x - a, y, z) on the edge where the isosurface intersects, M (x,y+b,z) , M (x,y-b,z) represents the gray values of (x, y + b, z) and (x, y - b, z), M (x,y,z+c) , M (x,y,z-c) represents the gray values of (x, y, z + c) and (x, y, z - c), and a, b, c respectively represent the spacings of two isosurface vertices in the X, Y, and Z axis directions; V x , V y , V z respectively represent the normal vectors of the triangular patch in the X, Y, and Z axis directions, G x1 , G y1 , G z1 are respectively the gradient components of point P1(x1, y1, z1) in the X, Y, and Z axis directions, G x2 , G y2 , Cz2 They are the gradient components of point P2(x2, y2, z2) in the X, Y, and Z axis directions respectively;

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

[0061] Define the n - order triangular domain Bezier surface as a triangular array composed of (n + 1)(n + 2) / 2 control vertices p i,j,k (i, j, k≥0, i + j + k = n):

[0062]

[0063] where P(u, v, w) is an arbitrary point on the triangular domain Bezier surface, u, v, w∈[0, 1] represent the barycentric coordinates inside the triangle, i, j, k represent the serial numbers of any three different control vertices, and B i,j,k n (u, v, w) represents the n - order Bezier basis function.

[0064] Calculate the barycentric coordinates inside the triangle through the following formula:

[0065]

[0066] where area represents the area of the triangle;

[0067] Convert the n - order control vertex P i,j,k to the (n - 1) - order control vertex Q i,j,k = uP i+1,j,k + vP i,j+1,k + wP i,j,k+1 , and recursively in turn. The finally remaining point is the point on the triangular domain Bezier surface, and thus reconstruct the root canal solid model.

[0068] Furthermore, perform coordinate alignment processing on the root canal solid model in 3D space, and use the approximate center line fitting the root canal to measure the morphological characteristics of the dental root canal, including:

[0069] Express the point cloud data x i =(x i , y i , z i )(i = 1, 2...n) of a single tooth as matrix X:

[0070]

[0071] where (x i , y i , z i) represents the coordinates of the i-th point cloud, and n represents the data volume of the point cloud;

[0072] The centering process of the point cloud data of a single tooth can be obtained through the following formula to get the matrix

[0073]

[0074] Among them,

[0075] Calculate through the following formula The covariance matrix of:

[0076]

[0077] In the formula, V represents The covariance matrix of, 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 dimensionally reduced, and D represents the covariance matrix of the white data;

[0078] If the eigenvalues of V are λ1, λ2,..., λ n , and the corresponding eigenvectors are p1, p2,..., p n , then there is:

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

[0080] It can be seen that the eigenvector corresponds to the rotation matrix R, representing the axis direction of each component; the eigenvalue corresponds to the stretching matrix S, representing the variance of the data in the axis direction of the corresponding component.

[0081] Arrange the eigenvectors in descending order of the corresponding eigenvalues from top to bottom to form a matrix and take the first m principal components to form the matrix W m , and determine that the Z-axis vector of the data principal component after the original point cloud is reduced to m dimensions is:

[0082]

[0083] Let the Z-axis vector of the principal component obtained according to formula (18) be k = (z x , z y , z z ), take i = (1, -z x / z y , 0) as the X-axis vector, and obtain the Y-axis vector by cross-multiplying the X-axis and the Z-axis to determine the principal component direction of the three-dimensional point cloud data of the tooth.

[0084] In a coordinate system aligned with the principal component direction, a method based on centerline approximation is used to measure the root canal length, the direction of the root canal orifice, the root canal curvature, and the pulp depth.

[0085] According to one aspect of the embodiments of the present application, a dental root canal feature morphology extraction device based on CBCT data is provided, including:

[0086] A three-dimensional reconstruction module, configured to perform filtering and denoising processing on the CBCT image and pulp segmentation processing, and perform surface rendering three-dimensional reconstruction on the dental root canal to generate a three-dimensional model;

[0087] A surface reconstruction module, 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;

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

[0089] According to one aspect of the embodiments of the present application, an electronic device is provided, including: a controller; a memory for storing one or more programs, when the one or more programs are executed by the controller, enabling the controller to implement the above-mentioned method for extracting the dental root canal feature morphology based on CBCT data.

[0090] According to one aspect of the embodiments of the present application, a computer-readable storage medium is further provided, on which computer-readable instructions are stored, when the computer-readable instructions are executed by a processor of a computer, enabling the computer to execute the above-mentioned method for extracting the dental root canal feature morphology based on CBCT data.

[0091] According to one aspect of the embodiments of the present application, a computer program product or a computer program is further provided, the computer program product or the computer program includes computer instructions, and the computer instructions are stored in a computer-readable storage medium. A processor of a computer device reads the computer instructions from the computer-readable storage medium, and the processor executes the computer instructions, enabling the computer device to execute the above-mentioned method for extracting the dental root canal feature morphology based on CBCT data.

[0092] In the technical solution provided by the embodiments of the present application, there are at least the following advantages:

[0093] This application is based on the CBCT tomographic images obtained by scanning a patient's teeth, and quickly reconstructs the corresponding three-dimensional model, which contains a large amount of morphological and positional information of the dental body, dental pulp, and root canal. Since the pulp chamber needs to be opened on the tooth crown surface before root canal preparation, obtaining relevant data such as the depth of the dental pulp and the size of the root canal orifice from the three-dimensional model of the tooth can help doctors determine the pulp chamber opening plan, which has practical significance and application value for the current root canal preparation surgery.

[0094] It should be understood that the above general description and the following detailed description are only exemplary and explanatory, and cannot limit this application. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] The accompanying drawings herein are incorporated into the specification and form a part of the specification, showing embodiments consistent with this application, and are used together with the specification to explain the principles of this application. Obviously, the accompanying drawings in the following description are only some embodiments of this application. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts. In the drawings:

[0096] Figure 1 is the overall flowchart of a method for extracting the characteristic morphology of tooth root canals based on CBCT data shown in an exemplary embodiment of this application;

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

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

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

[0100] Figure 5 is the overall flowchart of dental pulp segmentation shown in an exemplary embodiment of this application;

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

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

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

[0104] Figure 9It is a schematic diagram of the control network of a third-order triangular Bézier surface shown in an exemplary embodiment of the present application;

[0105] Figure 10 It is a schematic diagram of the recursive calculation process shown in an exemplary embodiment of the present application;

[0106] Figure 11 It is a flowchart of the measurement of the root canal morphology characteristics shown in an exemplary embodiment of the present application;

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

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

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

[0110] Figure 15 It is a structural diagram of a device for extracting the characteristic morphology of a root canal based on CBCT data shown in an exemplary embodiment of the present application. Detailed implementation manners

[0111] Here, the exemplary embodiments will be described in detail, and the examples are shown in the drawings. When the following description refers to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The implementation manners described in the following exemplary embodiments do not represent all implementation manners consistent with the present application. On the contrary, they are only examples of devices and methods consistent with some aspects of the present application as detailed in the appended claims.

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

[0113] The flowcharts shown in the drawings are only exemplary descriptions and do not necessarily include all contents and operations / steps, nor do they necessarily need to be executed in the described order. For example, some operations / steps can be decomposed, and some operations / steps can be combined or partially combined. Therefore, the actual execution order may change according to the actual situation.

[0114] As used in this application, "a plurality of" means two or more. " / or" describes the relationship between associated objects, indicating that there can be three relationships. For example, A / or B can represent: A exists alone, A and B exist simultaneously, and B exists alone. The character " / " generally indicates that the associated objects before and after are in an "or" relationship.

[0115] Please refer to Figure 1 , Figure 1 which is the overall flowchart of a method for extracting the characteristic morphology of dental root canals based on CBCT data shown in an exemplary embodiment of this application. One aspect of the embodiments of this application provides a method for extracting the characteristic morphology of dental root canals based on CBCT data. As Figure 1 shown, this method includes steps S100 to S300. The details are introduced as follows.

[0116] S100, perform filtering denoising and pulp segmentation processing on the CBCT image, and perform surface rendering 3D reconstruction on the tooth root canal to generate a 3D model.

[0117] In this embodiment, in order to reduce the interference of irrelevant tissues such as gums and muscles in the image, it is necessary to perform preprocessing operations of filtering denoising on the CBCT data before reconstructing the tooth surface. Then, perform surface rendering reconstruction on the preprocessed CBCT image to generate the 3D point cloud data of the target tooth.

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

[0119] S110, based on the obtained CBCT image, use the neighborhood edge model to improve the denoising effect of median filtering, and calculate the optimal threshold for segmenting the pulp based on the intercept histogram of reciprocal cross entropy to obtain the preprocessed CBCT sequence image.

[0120] In this embodiment, step S110 can be implemented through two steps. One step is CBCT image preprocessing, that is, based on the obtained CBCT image, use the neighborhood edge model to improve the denoising effect of median filtering; the other step is pulp segmentation, that is, calculate the optimal threshold for segmenting the pulp based on the intercept histogram of reciprocal cross entropy to obtain the preprocessed CBCT sequence image.

[0121] Specifically, when performing CBCT image preprocessing, first judge and retain the image edge. For example, in a 3×3 neighborhood, for the central pixel point M of this neighborhood c establish 8 edge models as Figure 3 shown, and calculate the pixel difference between this pixel point and its 8 neighboring points above, below, left, right, and at a 45° angle:

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

[0123] Wherein, M i represents the i-th neighboring point, i is a positive integer not greater than 8, and 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 value to 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] Wherein, 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, then 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 med, maximum value F max And repeat the first step. If the neighborhood area does not meet the second condition, the gray value of the central pixel point is not replaced; where 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 > 0, then the gray value of the central pixel point is not replaced. If not satisfied, the gray value of the central pixel point is replaced with the median value F med .

[0131] The comparison of the effects of using mean filtering and edge-preserving median filtering is as Figure 4 shown in (a) and (b) below.

[0132] The overall flowchart of pulp segmentation is as Figure 5 shown. When performing pulp segmentation, for a certain layer of image in CBCT, first use gamma transformation to enhance the image contrast, and then calculate the neighborhood average image and the gradient image respectively. Combine these two images with the gray image and add them pixel by pixel to form a mixed image. Based on the mixed image, its histogram can be used to find the optimal threshold. Then compare each pixel value in the mixed image with the optimal threshold and convert it into a binary image. Finally, clear all unnecessary connected tissues on the image boundary and only retain the pulp area.

[0133] Generate a three-dimensional space with gray information, neighborhood average gray information, and gradient compensation, as Figure 6 shown. The entire space is divided into a target region O and a background region B by a plane α.

[0134] Among them, the gray value of the point (x, y) and its neighborhood average gray value are respectively denoted as g(x, y) and f(x, y), 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 is the gradient compensation when the threshold is T. The gradient compensation for each threshold is a corresponding fixed value. Therefore, the optimal segmentation plane α is only determined by the threshold T, and the value of T is determined by the straight line k in the two-dimensional plane. Each line k intersects with the dotted line v, so the number of these lines is 2L - 1, that is, 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 respectively denoted as GH (x, y) and G V (x, y). The final gradient G(x, y) is a combination of the two directional gradients.

[0136]

[0137] In the formula, G H (x, y) and G V (x, y) represent the gradient values of the point (x, y) in the given image in the horizontal and vertical directions, x represents the abscissa value of the point in the given image, y represents the ordinate value of the point in the given image, G(x, y) represents the final gradient, g(x, y + 1) represents the gray value of the point (x, y + 1), g(x, y - 1) represents the gray value of the point (x, y - 1), g(x + 1, y) represents the gray value of the point (x + 1, y), and g(x - 1, y) represents the gray value of the point (x + 1, y);

[0138] Since the thresholds are the same and both the gray axis and the neighborhood average gray axis are T, adding the gray image, the neighborhood average image, and the gradient image generates a mixed image, expressed 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 the mixed image, the optimal threshold can be found by calculating the corresponding histogram. In the histogram, the mixed image with the threshold of T is divided into the target part Ω O ∈{(x, y)|F(x, y) = 0, 1, …, T} and the background part Ω B ∈{(x, y)|F(x, y) = T + 1, T + 2, …, 2T - 2} into two parts. Pixel values less than T belong to Ω O , and pixel values greater than T belong to Ω B .

[0142] Then, calculate the reciprocal cross-entropy between the target region and the background region through the following formula:

[0143]

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

[0145] The reciprocal cross-entropy E(O, B) can reflect the deviation before and after image segmentation. Therefore, when E(O, B) reaches the minimum value, there is:

[0146]

[0147] wherein, is the optimal threshold for dividing the dental pulp, and arg min is the operation of taking the index of the minimum element of the array.

[0148] S120. Draw the surface model of the root canal by surface rendering for the preprocessed CBCT sequence images.

[0149] The purpose of step S120 is to achieve 3D reconstruction, so as to obtain the surface model of the root canal, that is, the 3D dental model.

[0150] Exemplarily, uniform sampling is performed in the x, y, and z directions of the upper and lower jaw regions of the patient in 3D space, and the sampling intervals are Δx, Δy, and Δz. Then the volume data can be represented by the ternary function p(i, j, k). The cubic region formed by eight adjacent sampling points is a voxel, as Figure 7 shown.

[0151] The essence of the MC algorithm is to extract the isosurface for each voxel. The isosurface refers to the surface composed of points with the same attributes in 3D space. Different isosurface thresholds need to be set for reconstructing different tissues, and the threshold for teeth is set to about 750. Comparing the vertices of each voxel with the set isosurface threshold, these vertices all have two different states. When the states of the two vertices on one side of the voxel are different, this edge will definitely intersect with the isosurface. Due to the complementary symmetry and rotational symmetry of space, the configuration of the intersection of the voxel and the isosurface can be simplified to 15 types.

[0152] Traverse the entire volume data to find the voxels containing the isosurface, and use the method of linear interpolation to determine the intersection position P(x, y, z) of the isosurface and each edge of the voxel:

[0153]

[0154] wherein, T is the isosurface threshold, M1 and M2 are the gray values of the first vertex and the second vertex 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 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 calculated intersection points are the extracted isosurfaces. To render a better 3D effect for the isosurface triangular network, it is also necessary to select an appropriate lighting model according to the normal vectors of each triangular patch. First, calculate the gradients in the three-dimensional directions at the vertices of the isosurface:

[0156]

[0157] In the formula, G x , G y , G z respectively represent the gradients at the vertices of the isosurface in the X, Y, and Z axis directions,, M (x+a,y,z) , M (x-a,y,z) represent the gray values of two vertices (x + z, y, z) and (x - a, y, z) on the edge where the isosurface intersects, M (x,y+b,z) , M (x,y-b,z) represent the gray values of (x, y + b, z) and (x, y - b, z), M (x,y,z+c) , M (x,y,z-c) represent the gray values of (x, y, z + c) and (x, y, z - c), and a, b, c respectively represent the spacings between two isosurface vertices in the X, Y, and Z axis directions; V x , V y , V z are the normal vectors of the triangular patch in the X, Y, and Z axis directions respectively, G x1 , G y1 , G z1 are the gradient components of point P1(x1, y1, z1) in the X, Y, and Z axis directions respectively, G x2 , G y2 , G z2 are the gradient components of point P2(x2, y2, z2) in the X, Y, and Z axis directions respectively;

[0158] S200. Based on the regular point cloud data in the 3D model, the point-plane method is used to implement surface reconstruction to 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 patches on the root canal surface to achieve surface reconstruction of the root canal model.

[0160] Exemplarily, the triangular domain Bezier surface is composed of control vertices connected by elementary functions, that is, the nth-order triangular domain Bezier surface is composed of (n + 1)(n + 2) / 2 control vertices p i,j,k(where \(i, j, k \geq 0\) and \(i + j + k = n\)) defines a triangular array:

[0161]

[0162] In the formula, \(P(u, v, w)\) is an arbitrary point on the triangular domain Bezier surface, \(u, v, w \in [0, 1]\) represent the barycentric coordinates within the triangle, \(i, j, k\) represent any three different control vertices, and \(B\) i,j,k n (u, v, w) represents the \(n\)-th order Bezier basis function.

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

[0164] The barycentric coordinates within 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] During the calculation of the third-order triangular domain Bezier surface, 10 control vertices will be used, labeled as \(p\) i,j,k (where \(i, j, k \geq 0\) and \(i + j + k = 10\)). These 10 control vertices are connected by straight lines following the subscript order to generate a surface control network composed of triangles. The control vertices of the network correspond one by one to the nodes of the triangular domain as Figure 9 shown.

[0169] The \(n\)-th order control vertex \(P\) i,j,k will be converted into the \((n - 1)\)-th order control vertex \(Q\) i,j,k = \(uP\) i+1,j,k + \(vP\) i,j+1,k + \(wP\) i,j,j+1 . Recursively in this way, the final remaining point is the point on the triangular domain Bezier surface. The process of recursive calculation is as Figure 10 shown.

[0170] S300. Align the coordinate of the root canal solid model in the three-dimensional space, and measure the morphological characteristics of the dental root canal by using the approximate centerline that fits the root canal; among them, the morphological characteristics of the dental root canal include at least one of the root canal length, the direction of the root canal orifice, the root canal curvature, and the pulp depth.

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

[0172] Exemplarily, the measurement process of the root canal morphological characteristics is as Figure 11 shown, including the following steps S301 and S302.

[0173] Step S301, coordinate alignment.

[0174] The point cloud data x i =(x i , y i , z i )(i = 1, 2... n) of a single tooth is represented as matrix X:

[0175]

[0176] The matrix can be obtained by performing a centering process on the point cloud data of a single tooth through the following formula:

[0177]

[0178] where,

[0179] The covariance matrix of is calculated through the following formula:

[0180]

[0181] In the formula, V represents the covariance matrix of , 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 dimension-reduced and transformed, and D represents the covariance matrix of the white data;

[0182] If the eigenvalues of V are λ1, λ2,..., λ n , and the corresponding eigenvectors are p1, p2,..., p n , then there is:

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

[0184] It can be seen that the eigenvector corresponds to the rotation matrix R, representing the axis direction of each component; the eigenvalue corresponds to the stretching matrix S, representing the variance of the data in the axis direction of the corresponding component.

[0185] Arrange the eigenvectors in descending order of the corresponding eigenvalues from top to bottom to form a matrix, and take the first m principal components to form a matrix W m , and determine that the data after reducing the original point cloud to m dimensions is:

[0186]

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

[0188] Suppose the Z-axis vector of the principal component obtained according to formula (18) is k=(z x , z y , z z ), and take i=(1, -z x / z y , 0) as the X-axis vector, and obtain the Y-axis vector by cross-multiplying the X-axis and the Z-axis to determine the principal component direction of the three-dimensional tooth point cloud data.

[0189] Step S302, parameter calculation.

[0190] The centerline inside the tooth root canal is a space curve, and the projections of this centerline in the buccolingual and mesiodistal directions of the root canal are two plane curves. Therefore, first section the constructed root canal solid model on the two planes of the buccolingual and mesiodistal directions respectively, as Figure 12 shown, and two groups of cross-sectional curves of the root canal contour can be obtained.

[0191] Perform spline fitting on the horizontal midpoints of these two groups of cross-sectional curves respectively, and the projection curves of the root canal centerline in the buccolingual and mesiodistal directions can be obtained. Here, the XOZ plane represents the buccolingual direction, and the YOZ plane represents the mesiodistal direction, as Figure 13 shown.

[0192] Select the intersection of the two surfaces where these two projection curves are located and perpendicular to their corresponding buccolingual plane or mesiodistal plane, as Figure 14 shown, and the space curve formed after the intersection is the centerline of the tooth root canal.

[0193] From this, the length of the entire root canal and the opening direction of the root canal orifice can be calculated according to the approximate root canal center curve. The included angle between the tangent direction of a point at the root canal orifice of the centerline and the principal component direction of the tooth body is the opening direction of the root canal orifice.

[0194] Please refer to Figure 15 , Figure 15 which is the structural diagram of a device for extracting the characteristic morphology of tooth root canals based on CBCT data shown in an exemplary embodiment of the present application. On the other hand, the embodiment of the present application also provides a device for extracting the characteristic morphology of tooth root canals based on CBCT data. The device includes:

[0195] The three-dimensional reconstruction module 1501 is configured to perform filtering and denoising on CBCT images and pulp segmentation processing, and perform surface rendering three-dimensional reconstruction on tooth root canals 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 a three-dimensional space, and measure the morphological characteristics of the tooth root canal by using an approximate center line that fits the root canal; wherein, the morphological characteristics of the tooth root canal include at least one of root canal length, root canal orifice direction, root canal curvature, and pulp depth.

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

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

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

[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 adjacent points through the following formula:

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

[0204] In the formula, M i represents the i-th adjacent point, i is a positive integer not greater than 8, and the 8 adjacent 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 adjacent point;

[0205] Set the threshold size to T1, and record the number of adjacent points with gray levels similar to the central pixel point as n. If the pixel difference d between the central pixel point and the adjacent points is <T1, it is determined that the adjacent point has a gray level similar to the central pixel point; let n=n + 1. In the case of traversing all adjacent points, if n min ≤n≤n max , n minis the minimum number of neighboring points, n max is the maximum number of neighboring points, then determine that the number of neighboring points conforms to the established edge model, and determine that the central pixel point is an edge point and keep it directly;

[0206] For non-edge points, replace the non-edge points with the median of all pixel points in the neighborhood. Calculate the median of all pixel points through 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 is the median 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, M (x+m,y+n) is the pixel value of the pixel point (x+m,y+n);

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

[0210] 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;

[0211] The second step includes: expanding the neighborhood area by one circle. If the neighborhood area meets the second condition, sort all the pixel values in the neighborhood again, update the minimum value F min of all pixel points, the median value F med , and the maximum value F max and repeat the first step. If the neighborhood area does not meet the second condition, do not replace the gray value of the central pixel point; where 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 > 0, then do not replace the gray value of the central pixel point. If not satisfied, replace the gray value of the central pixel point 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 image, use gamma transformation to enhance the image contrast, and calculate the neighborhood average image and the gradient image respectively;

[0215] Combine the neighborhood average image and the gradient image with the grayscale image and add them pixel by pixel to form a mixed image;

[0216] Calculate the optimal threshold for segmenting the dental pulp using the intercept histogram of the mixed image;

[0217] Compare each pixel value in the mixed image with the optimal threshold and convert it into a binary image;

[0218] Clear all unnecessary connected tissues on the image boundary of the binary image, and only retain the dental pulp area to obtain the preprocessed CBCT sequence image;

[0219] Wherein, calculating the optimal threshold for segmenting the dental pulp using the intercept histogram of the mixed image includes:

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

[0221] Based on the three-dimensional space, calculate the final gradient through the following formula:

[0222]

[0223] In the formula, G H (x,y) and G V (x,y) represent the gradient values of the point (x,y) in the given image in the horizontal and vertical directions, x represents the abscissa value of the point in the given image, y represents the ordinate value of the point in the given image, G(x,y) represents the final gradient, g(x,y + 1) represents the grayscale value of the point (x,y + 1), g(x,y - 1) represents the grayscale value of the point (x,y - 1), g(x + 1,y) represents the grayscale value of the point (x + 1,y), and g(x - 1,y) represents the grayscale value of the point (x + 1,y);

[0224] Represent the mixed image 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, g(x,y) and f(x,y) represent the grayscale value of the point (x,y) and the neighborhood average grayscale value of the point (x,y);

[0227] Based on the segmentation threshold T, the mixed image is divided into the target part Ω O ∈{(x,y)|F(x,y) = 0, 1, …, T} and the 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] Calculate the reciprocal cross - entropy of the target region and the background region through the following formula:

[0229]

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

[0231] Based on the calculated reciprocal cross - entropy E(O,B), determine the optimal threshold for pulp segmentation through the following formula:

[0232]

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

[0234] In some embodiments, the 3D reconstruction module is further configured to:

[0235] When uniformly sampling in the x, y, and z directions of the upper and lower jaw regions in 3D space with sampling intervals of Δx, Δy, and Δz, represent the volume data through a ternary function; where a cube region composed of eight adjacent sampling points is a voxel;

[0236] Extract the isosurface for each voxel; wherein, the isosurface is a surface composed of points with the same attributes in 3D space;

[0237] Traverse the entire volume data, find the voxels containing the isosurface, and determine the intersection position P(x,y,z) of the isosurface and each edge of the voxel through the following formula:

[0238]

[0239] Wherein, T is the isosurface threshold, M1 and M2 are the gray values of the first vertex and the second vertex 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 side of the voxel;

[0240] Taking the triangular patch formed by the intersection points P(x, y, z) of the isosurface and each side of the voxel, 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 the isosurface, and combining each isosurface to form an isosurface triangular network;

[0241] In order to render a better three-dimensional effect of the isosurface triangular network, it is also necessary to select a lighting model according to the normal vectors of each triangular patch to render the isosurface triangular network to obtain the surface model of the root canal; wherein, the calculation formula for the normal vectors of each triangular patch is:

[0242]

[0243] Wherein, G x , G y , G z respectively represent the gradients in the X, Y, and Z axis directions at the isosurface vertex, M (x+a,y,z) , M (x-a,y,z) represent the gray values of two points (x + a, y, z) and (x - a, y, z) on the edge where the isosurface intersects, M (x,y+b,z) , M (x,y-b,z) represent the gray values of (x, y + b, z) and (x, y - b, z), M (x,y,z+c) , M (x,y,z-c) represent the gray values of (x, y, z + c) and (x, y, z - c), and a, b, c respectively represent the spacings of two isosurface vertices in the X, Y, and Z axis directions; V x , V y , V z respectively represent the normal vectors of the triangular patch in the X, Y, and Z axis directions, G x1 , G y1 , G z1 are respectively the gradient components of the point P1(x1, y1, z1) in the X, Y, and Z axis directions, G x2 , G y2 , G z2 are respectively the gradient components of the point P2(x2, y2, z2) in the X, Y, and Z axis directions;

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

[0245] The $n$-th order triangular domain Bézier surface is defined by a triangular array consisting of $\frac{(n + 1)(n + 2)}{2}$ control vertices $p$ i,j,k (where $i,j,k\geq0$ and $i + j + k = n$):

[0246]

[0247] where $P(u, v, w)$ is an arbitrary point on the triangular domain Bézier surface, $u, v, w\in[0, 1]$ represent the barycentric coordinates within the triangle, $i, j, k$ represent any three different control vertices, and $B$ i,j,k n (u, v, w) represents the $n$-th order Bézier basis function.

[0248] The barycentric coordinates within the triangle are calculated by the following formula:

[0249]

[0250] where area represents the area of the triangle;

[0251] The $n$-th order control vertex $P$ i,j,k is converted to the $(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 . Recursively, the final remaining point is the point on the triangular domain Bézier surface, and the root canal solid model is reconstructed accordingly.

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

[0253] The point cloud data $x$ i $=(x$ i , $y$ i , $z$ i ) ($i = 1, 2... n$) of a single tooth is represented as a matrix $X$:

[0254]

[0255] The point cloud data of a single tooth can be centered by the following formula to obtain the matrix

[0256]

[0257] where

[0258] The covariance matrix of is calculated by the following formula:

[0259]

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

[0261] If the eigenvalues of V are λ1, λ2,..., λ n , and the corresponding eigenvectors are p1, p2,..., p n , then there is:

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

[0263] It can be seen that the eigenvector corresponds to the rotation matrix R, representing the axis direction of each component; the eigenvalue corresponds to the stretching matrix S, representing the variance of the data in the axis direction of the corresponding component.

[0264] Arrange the eigenvectors in descending order of the corresponding eigenvalues from top to bottom to form a matrix, and take the first m principal components to form a matrix W m , and determine that the Z-axis vector of the principal component of the data after the original point cloud is reduced to m dimensions is:

[0265]

[0266] Suppose the Z-axis vector of the principal component obtained according to Equation (18) is k = (z x , z y , z z ), and use i = (1, -z x / z y , 0) as the X-axis vector, and obtain the Y-axis vector by cross-multiplying the X-axis and the Z-axis to determine the principal component direction of the three-dimensional tooth point cloud data.

[0267] Measure the root canal length, root canal orifice direction, root canal curvature, and pulp depth by the method based on the centerline approximation in the coordinate system aligned with the principal component direction.

[0268] It should be noted that the device for extracting the characteristic morphology of the tooth root canal based on CBCT data provided in the above embodiment and the method for extracting the characteristic morphology of the tooth root canal based on CBCT data provided in the foregoing embodiment belong to the same concept. The specific manner of performing the operations in the steps has been described in detail in the method embodiment and will not be repeated here.

[0269] Another aspect of the embodiment of the present application further provides an electronic device, including: a controller; a memory for storing one or more programs, and when the one or more programs are executed by the controller, to execute the method for extracting the characteristic morphology of the tooth root canal based on CBCT data in each of the above embodiments.

[0270] In particular, according to an embodiment of the present application, the processes described above with reference to the flowchart can be implemented as a computer software program. For example, an embodiment of the present application includes a computer program product that includes a computer program carried on a computer-readable medium, and the computer program includes a computer program for executing the method shown in the flowchart. In such an embodiment, the computer program can be downloaded and installed from the network through the communication part, and / or installed from a removable medium. When the computer program is executed by a central processing unit (CPU) 701, various functions defined in the system of the present application are executed.

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

[0272] The flowcharts and block diagrams in the accompanying drawings illustrate the possible architectures, functions, and operations of systems, methods, and computer program products according to various embodiments of the present application. Among them, each block in the flowchart or block diagram may represent a module, a program segment, or a part of code, and the above-mentioned module, program segment, or part of code contains one or more executable instructions for implementing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the blocks may occur in an order different from that marked in the accompanying drawings. For example, two consecutive blocks shown may actually be executed substantially in parallel, and they may sometimes be executed in the reverse order, depending on the functions involved. It should also be noted that each block in the block diagram or flowchart, as well as the combination of blocks in the block diagram or flowchart, can be implemented by a dedicated hardware-based system for performing the specified functions or operations, or can be implemented by a combination of dedicated hardware and computer instructions.

[0273] The modules / units involved in the embodiments described in the present application can be implemented in software or in hardware, and the described units can also be provided in the processor. Among them, the names of these modules / units do not constitute a limitation to the modules / units themselves in some cases.

[0274] Another aspect of the present application also provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, it implements the method for extracting the characteristic morphology of dental root canals based on CBCT data as described above. The computer-readable storage medium can be included in the electronic device described in the above embodiments, or can exist separately and not be assembled into the electronic device.

[0275] Another aspect of the embodiments of the present application also provides a computer program product or a computer program. The computer program product or the computer program includes computer instructions, and the computer instructions are stored in a computer-readable storage medium. The processor of the computer device reads the computer instructions from the computer-readable storage medium, and the processor executes the computer instructions, so that the computer device executes the method for extracting the characteristic morphology of dental root canals based on CBCT data provided in the above various embodiments.

[0276] According to one aspect of the embodiments of the present application, there is also provided a computer system, including a Central Processing Unit (CPU), which can perform various appropriate actions and processes according to a program stored in a Read-Only Memory (ROM) or a program loaded from a storage section into a Random Access Memory (RAM), such as executing the methods in the above embodiments. In the RAM, various programs and data required for system operation are also stored. The CPU, ROM, and RAM are connected to each other via a bus. An Input / Output (I / O) interface is also connected to the bus.

[0277] The following components are connected to the I / O interface: an input section including a keyboard, a mouse, etc.; an output section including, for example, a Cathode Ray Tube (CRT), a Liquid Crystal Display (LCD), etc. and a speaker, etc.; a storage section including a hard disk, etc.; and a communication section including a network interface card such as a LAN (Local Area Network) card, a modem, etc. The communication section performs communication processing via a network such as the Internet. A drive is also connected to the I / O interface as required. A removable medium, such as a magnetic disk, an optical disk, a magneto-optical disk, a semiconductor memory, etc., is installed on the drive as required so that a computer program read from it can be installed into the storage section as required.

[0278] The above content is only a preferred exemplary embodiment of the present application and is not used to limit the implementation of the present application. Those of ordinary skill in the art can make corresponding adaptations or modifications very conveniently according to the main concept and spirit of the present application. Therefore, the protection scope of the present application should be subject to the protection scope required by the claims.

Claims

1. A method for extracting characteristic morphology of root canals based on CBCT data, characterized in that: The method comprises: The CBCT images are filtered, de-noised and processed for dental pulp segmentation, and the tooth root canal is subjected to surface rendering and three-dimensional reconstruction to generate a three-dimensional model; Based on the regular point cloud data in the three-dimensional model, a point-to-surface method is used to realize surface reconstruction and generate a root canal solid model; The root canal entity model is aligned in three-dimensional space, and the root canal morphological characteristics are measured by fitting the approximate center line of the root canal; wherein the root canal morphological characteristics include at least one of the root canal length, root canal orifice direction, root canal curvature and pulp depth.

2. The method for extracting characteristic morphology of root canals based on CBCT data according to claim 1, characterized in that: The CBCT images are filtered, de-noised and segmented, and the root canals are reconstructed by surface rendering to generate a 3D model, including: Based on the acquired CBCT images, the neighborhood edge model is used to improve the denoising effect of the median filter, and the optimal threshold for segmenting the pulp is calculated based on the intercept histogram of the inverse cross entropy to obtain the preprocessed CBCT sequence images. The surface model of the root canal is drawn using the preprocessed CBCT sequence images through surface rendering.

3. The method for extracting characteristic morphology of root canals based on CBCT data according to claim 2, characterized in that: Based on the acquired CBCT images, the neighborhood edge model is used to improve the denoising effect of the median filter, including: In a neighborhood, eight edge models are established for the central pixel point of the neighborhood, and the pixel difference between the central pixel point and the neighboring points is calculated by the following formula: |M i -M c |=d(i=1,2,…,7,8) (1) Where M i represents the i-th neighboring point, i is a positive integer not greater than 8, and the 8 neighboring points are located above, below, left, right and in four 45° angles to the center pixel. c represents the central pixel, d represents the pixel difference between the central pixel and the i-th neighboring point; Set the threshold size to T1, and record the number of neighboring points with gray levels similar to the central pixel as n. If the pixel difference d between the central pixel and the neighboring points is < T1, then determine that the neighboring points have gray levels similar to the central pixel; let n = n + 1. After traversing all neighboring points, if n min ≤ n ≤ n max , n min is the minimum number of neighboring points, and n max is the maximum number of neighboring points, then determine that the number of neighboring 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 points are replaced with the median of all pixels in the 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 is the median of all pixels, and A is the center pixel ( x ,y) around the neighborhood, 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, M (x+m,y+n) is the pixel value of the pixel point (x+m,y+n); Sort all pixel values ​​in the neighborhood and get 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 >0, then execute the third step, if not, then execute the second step; The second step includes: expanding the neighborhood area, and if the neighborhood area meets the second condition, sorting all pixel values ​​in the neighborhood again, and updating the minimum value F of all pixel points min , median F med , maximum value F max Repeat the first step. If the neighborhood area does not meet the second condition, the gray value of the central pixel is not replaced; 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 >0, the grayscale value of the central pixel is not replaced; if not satisfied, the grayscale value of the central pixel is replaced with the median F med .

4. The method for extracting characteristic morphology of root canals based on CBCT data according to claim 2, characterized in that: The optimal threshold for segmenting the dental pulp is calculated based on the intercept histogram of the inverse cross entropy, and the preprocessed CBCT sequence images are obtained, including: Based on the denoised CBCT images, the image contrast is enhanced by gamma transformation, and the neighborhood average image and gradient image are calculated respectively. The neighborhood average image and the gradient image are combined with the grayscale image and added pixel by pixel to form a mixed image; The optimal threshold for segmenting the dental pulp was calculated using the intercept histogram of the mixed image; Compare each pixel value in the mixed image with the optimal threshold and convert it into a binary image; All unnecessary connected tissues on the image boundary of the binary image are removed, and only the pulp region is retained to obtain a preprocessed CBCT sequence image; The method of calculating the optimal threshold for segmenting the dental pulp by using the intercept histogram of the mixed image includes: The grayscale information, the neighborhood average grayscale information and the gradient compensation generate a three-dimensional space, wherein the three-dimensional space is divided into a target area O and a background area B by a plane; Based on the three-dimensional space, the final gradient is calculated by the following formula: In the formula, G H (x,y) and G V (x,y) represents the gradient value of the point (x,y) in the given image in the horizontal and vertical directions, x represents the abscissa value of the point in the given image, y represents the ordinate value of the point in the given image, G(x,y) represents the final gradient, g(x,y+1) represents the grayscale value of the point (x,y+1), g(x,y-1) represents the grayscale value of the point (x,y-1), g(x+1,y) represents the grayscale value of the point (x+1,y), and g(x-1,y) represents the grayscale value of the point (x+1,y); The mixed image is represented as: F(x,y)=f(x,y)+g(x,y)-G(x,y) (5) Where F(x,y) represents the mixed image, 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 mixed image is divided into the target part Ω O ∈{(x,y)|F(x,y)=0,1,…,T} and background part Ω B ∈{(x,y)|F(x,y)=T+1,T+2,…,2T-2}, the pixel values ​​less than the segmentation threshold T belong to the target part, and the pixel values ​​greater than the segmentation threshold T belong to the background part; The inverse cross entropy between the target area and the background area is calculated by the following formula: Where E(O,B) represents the inverse cross entropy; P(k) represents the prior probability that the average gray value of the neighborhood of the mixed image is k. represents the grayscale mean of the target area, Represents the grayscale mean of the background area; Based on the calculated cross entropy E(O,B), the optimal threshold for segmenting the pulp is determined by the following formula: In the formula, is the optimal threshold for segmenting the dental pulp, and arg min is the operation of taking the minimum element index of the array.

5. The method for extracting characteristic morphology of root canals based on CBCT data according to claim 4, characterized in that: The surface model of the root canal is drawn by surface rendering using the preprocessed CBCT sequence images, including: In the case where the x, y, and z directions of the upper and lower jaw regions are uniformly sampled in three-dimensional space and the sampling interval is Δx, Δy, Δz, the volume data is represented by a ternary function; a cubic area formed by eight adjacent sampling points is a voxel; Extracting an isosurface for each voxel; wherein the isosurface is a surface composed of points with the same attributes in three-dimensional space; Traverse the entire volume data to find the voxels containing the isosurface, and determine the intersection position P(x, y, z) of the isosurface and each edge of the voxel using the following formula: Where T is the isosurface threshold, M1, M2 are the grayscale 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 on the isosurface and each edge of the voxel; The intersection points P(x, y, z) on the isosurface and each edge of the voxel, and the triangular face formed by the first vertex P1(x1, y1, z1) and the second vertex P2(x2, y2, z2) on the edge where the voxel and the isosurface intersect are used as isosurfaces, and each isosurface is combined to form an isotriangulated network; In order to render the equivalent triangular network with a better three-dimensional effect, it is also necessary to select a lighting model according to the normal vector of each triangular facet, so as to render the equivalent triangular network to obtain a surface model of the root canal; wherein the calculation formula of the normal vector of each triangular facet is: In the formula, G x ,G y ,G z Respectively represent the gradient of the isosurface vertex in the X, Y, and Z axis directions, M (x+a,y,z) , M (x-a,y,z) Represents the grayscale values ​​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) Represents the grayscale value of (x,y+b,z) and (x,yb,z), M (x,y,z+c) , M (x,y,z-c) Represents the grayscale value of (x,y,z+c) and (x,y,zc), a, b, c represent the distance between the two isosurface vertices in the X, Y, Z axis directions respectively; V x ,V y ,V z The normal vectors of the triangle in the X, Y, and Z axis directions, G x1 ,G y1 ,G z1 are the gradient components of point P1 (x1, y1, z1) in the X, Y, and Z axis directions, G x2 ,G y2 ,G z2 They are the gradient components of point P2 (x2, y2, z2) in the X, Y, and Z axis directions respectively.

6. The method for extracting characteristic morphology of root canals based on CBCT data according to claim 1, characterized in that: Based on the regular point cloud data in the three-dimensional model, a point-to-surface method is used to realize surface reconstruction and generate a root canal solid model, including: The n-order triangular Bezier surface consists of (n+1)(n+2) / 2 control vertices p i,j,k The definition of the triangular array composed of (i, j, k ≥ 0, i + j + k = n) is: Where P(u,v,w) is any point on the Bezier surface of the triangle domain, u,v,w∈[0,1] represents the coordinates of the center of gravity in the triangle, i,j,k represents the serial numbers of any three different control vertices, and B i,j,k n (u,v,w) represents the n-th order Bezier basis function. The coordinates of the center of gravity in the triangle are calculated using the following formula: In the formula, area represents the area of ​​the triangle; Set the n-th order control vertex P i,j,k Convert to n-1 order control vertex Q i,j,k =uP i+1,j,k +vP i,j+1,k +wP i,j,k+1 , recursively, and finally the remaining points are points on the triangular domain Bezier surface, thereby reconstructing the root canal solid model.

7. The method for extracting characteristic morphology of root canals based on CBCT data according to claim 1, characterized in that: The root canal entity model is aligned with coordinates in three-dimensional space, and the root canal morphological characteristics are measured by fitting the approximate center line of the root canal, including: The point cloud data of a single tooth x i =(x i ,y i ,z i )(i=1,2...n) is represented as matrix X: Calculated by the following formula The covariance matrix of is: In the formula, V represents The covariance matrix of , S represents the point cloud stretching matrix, satisfying S = S T , R represents the point cloud rotation matrix, L represents the covariance of the new coordinate system after the point cloud is transformed into a dimensionality-reduced state, and D represents the covariance matrix of the white data; If the eigenvalues ​​of V are λ1,λ2,...,λ n , the corresponding eigenvectors are p1,p2,...,p n , then: V·p i =λ i ·V (17) It can be seen that the eigenvector corresponds to the rotation matrix R, which represents the coordinate axis direction of each component; the eigenvalue corresponds to the stretching matrix S, which represents the variance of the data in the coordinate axis direction of the corresponding component. Arrange the eigenvectors into a matrix from top to bottom according to the corresponding eigenvalue size and take the first m principal components to form the matrix W m , and determine the Z-axis vector of the principal component of the data after the original point cloud is reduced to m dimensions: Suppose the Z-axis vector of the principal component obtained according to formula (18) is k=(z x ,z y ,z z ), with i = (1, -z x / z y ,0) as the X-axis vector, and the Y-axis vector is obtained by multiplying the X-axis by the Z-axis to determine the principal component direction of the tooth 3D point cloud data. Root canal length, root canal orifice direction, root canal curvature, and pulp depth were measured using a centerline approximation method in a coordinate system aligned with the principal component directions.

8. A root canal feature morphology extraction device based on CBCT data, characterized in that: The device comprises: A three-dimensional reconstruction module is configured to perform filtering and 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; A surface reconstruction module is configured to realize surface reconstruction based on the regular point cloud data in the three-dimensional model by using a point-to-surface method to generate a 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 characteristics by fitting the approximate center line of the root canal; wherein the root canal morphological characteristics include at least one of the root canal length, root canal orifice direction, root canal curvature and pulp depth.

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

10. A computer-readable storage medium, characterized in that: Computer-readable instructions are stored thereon, and when the computer-readable instructions are executed by a processor of a computer, the computer is caused to execute the method for extracting characteristic morphology of root canals based on CBCT data according to any one of claims 1 to 7.

Images