CN110942433A - Skull CBCT image-based repair guide plate generation method - Google Patents
Skull CBCT image-based repair guide plate generation method Download PDFInfo
- Publication number
- CN110942433A CN110942433A CN201911148810.0A CN201911148810A CN110942433A CN 110942433 A CN110942433 A CN 110942433A CN 201911148810 A CN201911148810 A CN 201911148810A CN 110942433 A CN110942433 A CN 110942433A
- Authority
- CN
- China
- Prior art keywords
- skull
- guide plate
- curvature
- image
- points
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 65
- 210000003625 skull Anatomy 0.000 title claims abstract description 60
- 238000007408 cone-beam computed tomography Methods 0.000 title claims abstract description 37
- 230000008439 repair process Effects 0.000 title claims abstract description 16
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 31
- 238000000605 extraction Methods 0.000 claims abstract description 20
- 230000007547 defect Effects 0.000 claims abstract description 18
- 230000008569 process Effects 0.000 claims abstract description 15
- 239000013598 vector Substances 0.000 claims description 23
- 239000011159 matrix material Substances 0.000 claims description 7
- 210000003484 anatomy Anatomy 0.000 claims description 5
- 238000010606 normalization Methods 0.000 claims description 4
- 238000009499 grossing Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000003287 optical effect Effects 0.000 claims description 2
- 238000007781 pre-processing Methods 0.000 claims description 2
- 238000004458 analytical method Methods 0.000 claims 1
- 238000003745 diagnosis Methods 0.000 abstract description 3
- 238000004364 calculation method Methods 0.000 description 9
- 238000010586 diagram Methods 0.000 description 6
- 210000000988 bone and bone Anatomy 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000001356 surgical procedure Methods 0.000 description 5
- 238000013461 design Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 230000011218 segmentation Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000010146 3D printing Methods 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000036772 blood pressure Effects 0.000 description 1
- 230000036770 blood supply Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000002316 cosmetic surgery Methods 0.000 description 1
- 238000012938 design process Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000003708 edge detection Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000036244 malformation Effects 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 210000000214 mouth Anatomy 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000399 orthopedic effect Effects 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 239000011800 void material Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/77—Retouching; Inpainting; Scratch removal
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30008—Bone
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30016—Brain
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Image Processing (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The invention discloses a skull CBCT image-based repair guide plate generation method, which comprises the steps of detecting skull defect edge point cloud through an edge extraction method, calculating the curvature of a defect region in a mirror image, further executing a region growing algorithm through the curvature to form a first-stage repair guide plate, and finally optimizing the first-stage repair guide plate to improve the consistency of an image and an actual form. The skull repairing structure pre-judging method is simple in operation flow, and a doctor performs pre-judgment on the skull repairing structure in advance in the digital preoperative planning and diagnosis process, so that the preoperative planning and the intraoperative conditions are consistent, the planning accuracy of the operation is quantized, the consciousness of the digitization degree of the doctor is improved, and the user experience is enhanced.
Description
Technical Field
The invention relates to the technical field of skull surgeries, in particular to a repair guide plate generation method based on skull CBCT images.
Background
In the process of reshaping the skull by reconstructing the skull top, the most critical technology is to redesign the original skull, namely to cut, excise and split the skull, so as to achieve the effect of reshaping. The current skull top design judges the skull top malformation condition by three-dimensional CT, designs the cutting shape and the cutting route on a two-dimensional CT sheet or paper by the experience of a doctor, and then adjusts the cutting shape, the cutting route and the splicing mode on the spot according to the actual condition in the operation process.
However, the above conventional methods have many problems: 1. some skull top reconstruction does not necessarily need to cut and splice all the skull, and only part of the skull is cut, separated and spliced, so that the design of a scheme before the operation is difficult to make; 2. the skull top reconstruction relates to the positioning problem of a three-dimensional space, and the three-dimensional space relation among bone flaps in the skull top reconstruction process is difficult to completely consider on a two-dimensional plane; 3. due to the lack of effective means, the surgical scheme cannot be effectively transferred to the implementation process of the surgery, the cut bone flaps cannot be optimally utilized, and the accurate splicing relation between the bone flaps is difficult to obtain, so that the expected plastic result designed by the preoperative scheme is difficult to achieve; 4. because the surgical plan cannot be fully implemented in the surgical procedure, the protection of important positions and anatomical structures planned in the preoperative plan cannot be realized in the surgical procedure; 5. the skull bone cutting instrument is lack of effective auxiliary tools, and cannot be meticulous and accurate in the operation process, so that the skull bone can be excessively cut, bone flaps are greatly damaged, and the bone flaps are not easy to heal after the operation; 6. the operation time is long, the blood supply, the blood pressure and the metabolism of a patient are influenced, uncertain factors of the operation are increased, and the operation risk is increased.
At present, the digital biotechnology is rapidly developed and widely applied to surgical clinics. Among them, the techniques that can accurately transfer preoperative designs to intraoperative operations include surgical navigation and guide plates. The operation navigation needs to match and collate the information of CT, MRI and the like before or during the operation with the entity data of the patient through a large-scale and expensive navigator, so that the operation is complex, and the possibility of navigation failure caused by matching error exists; the surgical guide plate technology has been applied to orthopedics, plastic surgery and maxillofacial surgery to a certain extent, but most of the used patient data are from CT or three-dimensional scanning, the requirement on equipment is high, the guide plate design process is also complex, special engineering designers are needed, and the clinical popularization is not facilitated.
Therefore, in the prior art, reconstructed volume data of the skull is obtained after CBCT reconstruction, so that the morphological structure of the skull in a three-dimensional space is analyzed, relevant information is obtained, and a panoramic image of teeth can be generated by the reconstructed volume data, so that X-ray panoramic scanning can be performed on an object less, and the radiation dose is reduced. However, most products in the prior art cannot support the function of root checking based on oral CBCT images, and cannot achieve the positioning effect of digital diagnosis in 3D views; meanwhile, in the operation performed by scanning methods such as CT, the matching degree of images and entities is low, the requirement of the skull operation on the accuracy is very high, if the relative position in the skull cannot be judged during the operation, the operation cannot be performed, and further the condition that a doctor cannot pre-judge the skull in the digital preoperative planning process is caused, so that the risk coefficient of the operation and the failure rate of the operation are increased.
Disclosure of Invention
The invention aims to provide a guide plate repairing method based on a skull CBCT image, which improves the matching consistency of skull reconstruction data and patient entity data, improves the extraction precision, and improves the safety factor and the success rate of the operation.
The technical scheme adopted by the invention for solving the technical problems is as follows: a skull-based CBCT image restoration guide plate generation method is characterized by comprising the following steps:
s100, obtaining a three-dimensional skull model through optical scanning, and establishing a projection image and a CBCT file sequence;
s200, analyzing a CBCT file sequence through a Marching Cube and a region growing algorithm;
s300, detecting a skull defect edge point cloud by an OpenCV edge extraction method, and generating a mirror image by using symmetrical points in the point cloud;
s400, calculating the curvature of the defect area in the mirror image, wherein the curvature comprises an average curvature, a Gaussian curvature and a main curvature;
s500, executing a region growing algorithm based on the curvature in the S400 to form a first-stage repair guide plate;
s600, optimizing the edge of the first-stage repairing guide plate based on a Poisson equation hole-filling smoothing method to form a final repairing guide plate.
Preferably, the Marching Cube parsing process includes the following steps:
s201, reading the original data into a specified array after preprocessing;
s202, extracting a unit body from the grid data body to form a current unit body, and simultaneously acquiring all information of the unit body, wherein the information comprises a vertex function value of a boundary and a point cloud coordinate position of the unit body;
s203, comparing the function value of the vertex in the current unit body with the given isosurface value C to obtain a state table of the unit body;
s204, finding out the unit body edge intersected with the isosurface according to the current state table index of the unit body, and calculating the position coordinates of each intersection point by adopting a linear interpolation method;
s205, solving the normal vector of the vertex in the current unit body by using a center difference method, and obtaining the normal vector of each vertex of the triangular patch by using a linear interpolation method;
and S206, drawing the isosurface image according to the coordinates of the vertexes of the triangular patches and the vertex normal vectors.
Preferably, the executing of the region growing algorithm in S200 includes controlling the threshold value at 770-25584 intervals, and reconstructing a 3D model to view the main anatomical structure of CBCT.
Preferably, the OpenCV edge extraction method in S300 includes the following steps:
generating a group of normalized Gaussian kernels through a discretized Gaussian function, and then carrying out weighted summation on each point of a gray matrix in the image data based on the Gaussian kernels;
highlighting points with significant changes in the intensity values of the neighborhood of the gray points in the image data through an enhancement algorithm;
and eliminating the noise edge points by a thresholding method.
Preferably, the mirror image generated by the symmetry point in S300 includes an extraction of feature points, and an extraction formula thereof is as follows:
wherein p isi,pjAre two points of symmetry, hi,hjIs the result after normalization, Ci,jIs the vector obtained, C is extractedi,jMaximum value to obtain feature points.
Preferably, the mean curvature of the defect region in S400 is obtained by the following formula:
wherein, P and Q are curvature mean values of the characteristic points, C is a form distance, and T is TPS estimated transformation matrix weight.
Preferably, the S500 comprises a region growing algorithm, which comprises the following steps:
s501, selecting a seed point p (x)0，y0) The stack represents a seed area, and the seed point is pushed into the seed stack;
s502, stacking a first seed point pop in the seed stack, and traversing pixels in the neighborhood of the center 8 by taking the point as the center;
s503, judging whether the traversal pixel point is in the seed region, if not, judging whether the traversal pixel point meets the similarity of adjacent seed points, and if the pixel point (x, y) meets the similarity, pushing the pixel point (x, y) into a stack;
s504, repeating S502-S503 until the seed stack is empty.
Preferably, the step S600 includes generating an initialized hole filling grid according to the hole boundaries in the triangular grid, and then modifying the geometry of the triangular patch in the hole filling grid through normal estimation and Poisson equation, so that the geometry can be adapted and fused with the surrounding original grid, and the process includes the following steps:
s601, detecting hole boundaries and initializing hole filling grids;
s602, calculating an expected normal direction of the vertexes in the hole filling grid, constructing a Laplace equation delta f which is 0, and solving normal distribution of the vertexes inside the hole filling grid, wherein a Laplace operator is as follows:
wherein N is1(xi) Representing a vertex xi1 ring neighborhood point of αijAnd βijIs an edge eij2 corresponding diagonal corners;
s603, rotating the triangular patches in the hole filling grid based on the expected normal direction;
s604, adjusting the vertex position of the hole filling grid based on a Poisson equation, and calculating the gradient field of the three-dimensional grid, wherein the formula is as follows:
f is the position of the peak of the adjusted grid to be solved, and w is the gradient field of the torn grid;
the gradient operator is:
wherein the expression of the basis function gradient ▽ Φ i is:
it represents the vector rotated 90 degrees counterclockwise, AT represents the area of the triangle T;
divergence operator:
wherein, T1(xi) Representing a vertex xiAnd (3) the 1-ring neighborhood triangle, AT represents the area of the triangle T.
And S605, performing iterative region growing based on the curvature of the defect region obtained in the S400 to obtain a 3D point cloud patch.
The invention has the beneficial effects that: the invention forms a simulation image by extracting the skull morphology, helps doctors predict skull trends in and after the operation in advance, greatly improves the accuracy of the operation and reduces the operation risk. Meanwhile, the skull defect edge point cloud is detected by an edge extraction method, the curvature of a defect area in a mirror image is calculated, a region growing algorithm is executed through the curvature to form a first-stage repair guide plate, and finally, optimization processing is carried out on the first-stage repair guide plate to improve the consistency of the image and the actual form. The simulation image can be used for digital simulation preoperative planning, and can be used for preoperative simulation testing through a 3D printing generation entity, and the situations occurring in the operation can be displayed in advance in a simulation mode. The skull repairing structure pre-judging method is simple in operation flow, and a doctor performs pre-judgment on the skull repairing structure in advance in the digital preoperative planning and diagnosis process, so that the preoperative planning and the intraoperative conditions are consistent, the planning accuracy of the operation is quantized, the consciousness of the digitization degree of the doctor is improved, and the user experience is enhanced. In addition, through the test of the inventor, the processing method is also applicable to the extraction process of the fine anatomical structures in CBCT images and MRI images related to other medical departments.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the present invention will be further described with reference to the accompanying drawings and embodiments, wherein the drawings in the following description are only part of the embodiments of the present invention, and for those skilled in the art, other drawings can be obtained without inventive efforts according to the accompanying drawings:
FIG. 1 is a schematic view of a scanning model of step S100 in a skull-based CBCT image restoration guide generation method according to the present invention;
FIG. 2 is a schematic diagram of a skull with a defect area in step S100 of the method for generating a repair guide plate based on a CBCT image of the skull according to the present invention;
FIG. 3 is a schematic mirror image generated by using symmetric points in a point cloud in step S200 in the skull-based CBCT image repairing guide plate generating method of the present invention;
FIG. 4 is a schematic curved surface diagram of a defect area in a mirror image in a skull CBCT image-based repair guide plate generation method of the present invention;
FIG. 5 is a schematic diagram of a first-stage repairing guide plate in the repairing guide plate generating method based on the skull CBCT image of the invention;
FIG. 6 is a schematic diagram of a final repair guide plate in a mirror image in a repair guide plate generation method based on a skull CBCT image according to the present invention;
FIG. 7 is a schematic diagram of a point cloud patch in a skull-based CBCT image restoration guide plate generation method of the present invention;
FIG. 8 is a schematic diagram of a Gaussian-Boeing-internal-treatment result in the skull-based CBCT image restoration guide generation method of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the following will clearly and completely describe the technical solutions in the embodiments of the present invention, and it is obvious that the described embodiments are some embodiments of the present invention, but not all embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments of the present invention without inventive step, are within the scope of the present invention.
In embodiment 1, as shown in fig. 1, 2, 3, 4, 5 and 6, a method for generating a repair guide based on a CBCT image of a skull, the method comprising the steps of: firstly, analyzing a file sequence (dcm format) of the CBCT in the oral cavity, executing a region growing algorithm, and simultaneously controlling the CBCT image threshold value in 770-25584 intervals to further reconstruct a main anatomical structure (3D view) of the CBCT; detecting a skull defect edge point cloud by an OpenCV edge extraction method, and generating a mirror image by using symmetrical points in the point cloud; then calculating the curvature of the defect area in the mirror image, and executing an area growing algorithm through the curvature to form a first-stage repair guide plate; and finally, optimizing the edge of the first-stage repaired guide plate based on a Poisson equation hole filling smoothing method to form the final repaired guide plate.
Specifically, for the Marching Cube method and region growing algorithm of step S200, Marching Cube utilizes an iso-surface, i.e., a collection of all points in three-dimensional stereo space having the same value, whose collection region exhibits contours resembling those in a topographic map. Considering the basic assumption in Marching Cube that there is a Cube (hexahedron) in space and the data field along the edges of the hexahedron is continuous, that is, when two vertices of an edge of the Cube are respectively larger or smaller than the value of the isosurface, there is and only one point on the edge is the intersection of the edge with the isosurface. In the process of the space, namely, a large number of small cubes are used for segmenting the space, and then the plane inside the small cubes is used for approximately representing the current isosurface, so that the larger the number of the small cubes is, the better the approximation effect is, and a large number of calculations are carried out accordingly. The calculation formula comprises { (x, y, z) | f (x, y, z) ═ c }, wherein c is a constant, cubes (voxels) in the data field are processed one by one, the calculation process comprises the steps of separating the cubes intersected with the isosurface, calculating the intersection points of the isosurface and the edges of the cubes by adopting bilinear interpolation, and connecting the intersection points of the isosurface and the edges of the cubes vertically in a three-dimensional space according to the relative position of each vertex of the cubes and the isosurface to generate the isosurface, so that the isosurface is used as an approximate representation of the isosurface in the cubes.
Further, step S201 is executed, the original data is read into a designated array after being preprocessed; step S202 is executed to extract a unit from the mesh data volume to become a current unit, and acquire all information of the unit, where the information includes a function value of a vertex of the boundary and a point cloud coordinate position of the unit, and in this embodiment, 8 vertices are taken as an example; step S203 is executed, the function values of the 8 vertexes of the current unit body are compared with the given isosurface value C, and the state of the unit body is obtained; (edgeTable, triTable); step S204 is executed, the unit body edge intersected with the isosurface is found out according to the current state table index of the unit body, and the position coordinates of each intersection point are calculated by adopting a linear interpolation method; step S205 is executed, normal vectors of 8 vertexes of the current unit body are obtained by using a center difference method, and the normal direction of each vertex of the triangular patch is obtained by adopting a linear interpolation method; and step S206 is executed, and the isosurface image is drawn according to the coordinates of the vertexes of the triangular patches and the vertex normal vectors.
Specifically, in the center difference method calculation process, for example, if equal time steps are used, Δ t (i) ═ Δ t (Δ t is a constant), and u represents the displacement, then the center difference of the velocity and the acceleration is approximated by: u' (i) ═ u (i +1) -u (i-1) ]/(2 Δ t); u "(i) ═ u (i +1) -2u (i) + u (i-1) ]/(Δ t × Δ t).
For the region growing algorithm in step S200, a seed point (i.e., a feature point) in the partition region is obtained, pixel points having similar properties to the seed point are searched around the seed point and merged into the seed region, and then the merged pixels are continuously searched as new seed points until there are no similar pixel points around all pixels in the seed region, and the algorithm is ended. The basic algorithm flow is as follows: step 1, selecting a seed point p (x)0，y0) The stack represents a seed area, and the seed point is pushed into the seed stack; step 2, stacking a first seed point pop in the seed stack, and traversing pixels in the neighborhood of the center 8 by taking the point as the center; judging whether a traversal pixel point is already in the seed region, if not, judging whether the traversal pixel point meets the similarity of adjacent seed points, and if the pixel point (x, y) meets the similarity, pushing the pixel point (x, y) into a stack; and 4, repeating the steps 2 to 3 until the seed stack is empty. It can be seen that there are three elements that affect the region growing algorithm: selecting seed points; selecting a search path; and (5) judging the pixel similarity.
For the selection of the seed points, the region growing algorithm is a semi-interactive segmentation algorithm, and the seed points are required to be selected by a user or calculated by other algorithms; for the selection of the search path, the search path is to select adjacent pixels, taking a two-dimensional image as an example, and generally performing 8-neighborhood search, or 4-neighborhood search, taking a three-dimensional image as an example, and generally performing 26-neighborhood search or 6-neighborhood search; for the judgment of the pixel similarity, the similarity is judged by taking the degree of similarity of the pixel values as a judgment standard, for example, a certain gray scale range is set as a similarity standard, or a certain shape or property can be satisfied by calculation as a judgment standard.
Further, in the OpenCV edge extraction method in step S300, the extraction steps are as follows:
filtering:
the algorithms for edge detection are mainly based on the first and second derivatives of the image intensity, but the derivatives are usually very sensitive to noise, so filters must be used to improve the performance of the noise-related edge detector. For example, filtering is performed by using a gaussian filtering method, that is, a set of normalized gaussian kernels is generated by using a discretized gaussian function, and then each point of the image gray matrix is subjected to weighted summation based on the gaussian kernels.
Enhancing:
the basis of the enhanced edge is to determine the variation value of the neighborhood intensity of each point of the image. The enhancement algorithm can highlight points with significant changes in the intensity values of the image gray point neighborhood. In a particular programming implementation, this is determined by calculating the corresponding gradient magnitude.
And (3) detection:
the enhanced image has a large number of gradient values in its neighborhood, which are not the edge points to be found in a particular application, and therefore need to be rounded off. In actual engineering, edge points can be detected by a thresholding method to extract the edge points.
In this embodiment, edge extraction is performed by a method based on a Laplacian operator-Laplacian function, and the function call form is as follows:
void Laplacian(inputArray,outputArray,int ddepth,int ksize＝1,doublescale＝1,double delta＝0,int borderType＝BORDER_DEFAULT)。
further, the mirror image generated by the symmetric point in S300 includes extraction of feature points, and an extraction formula thereof is as follows:
wherein p isi,pjAre two points of symmetry, hi,hjIs the result after normalization, Ci,jIs the vector obtained, C is extractedi,jMaximum value to obtain feature points.
And for the selection of the symmetrical points, the two farthest eyepit characteristic points and the vertical normal line part of the nose bridge in the skull are mainly and automatically selected, the point cloud matrix is automatically copied at the same time, the triangulation is carried out to generate a mirror image, and the maximum value of Ci and j in an eyepit vector is obtained in a 3D space to obtain the characteristic points.
Further, the curvature in step S400 includes the following steps of mean curvature, gaussian curvature and principal curvature:
the principle of calculation of the mean curvature is as follows:
the average curvature is related to the unit normal vector of the surface, the selection of which affects the sign of the curvature. The sign of the curvature depends on the direction of the normal vector: the curvature is positive if the surface is "far" from the normal vector. The above definition holds for any defined surface in 3D space, which is mainly used to calculate the divergence of unit normal vectors.
For a surface defined as a function of two coordinates, e.g. z ═ S (x, y), expressed as (twice) the average curvature using a downward normal vector
The principle of gaussian curvature is calculated as follows:
the sum of the triangles of the triangle on the negative curvature surface is less than the sum of the triangles of the plane triangle.
The total curvature of the geodetic triangle (i.e. the triangle in the riemann sphere geometry) is equal to the difference of the sum and of its internal angles. The sum of the internal angles of the triangles on the positive curvature surface is greater than 180 degrees, and the sum of the internal angles of the triangles on the negative curvature surface is less than 180 degrees. On a zero curvature surface (e.g., a euclidean plane), the sum of its internal angles is equal to pi.
A more general result, as shown in figure 1, is the gaussian-blonde theorem. It can be seen that the sum of the triangles of the triangle on the negative curvature surface is smaller than the sum of the triangles of the plane triangle.
The principle of calculation of the principal curvatures is as follows:
for a given point P0(u0, v0) on the curved surface S r (u, v), the normal curvature kn is a function of the tangent direction du dv, and each critical value of the normal curvature is called the main curvature of the curved surface at this point; the corresponding direction is referred to as the main direction of the curved surface at this point. The following theorem was used for the calculation:
one point on the curved surface is represented by equation Pdu2+2Qdudv+Rdv2The two tangential directions perpendicular to each other defined as 0 are the conditions of ER-2FQ + GP 0
Wherein E, F, G is the first type basis quantity for a curved surface. The requirement for orthogonality between the two directions du dv and δ u δ v is Edu δ u + F (du δ v + dv δ u) + Gdv δ v ═ 0, i.e.:
writing a known quadratic equation as
Wherein the content of the first and second substances,
is two roots thereof, and both should satisfy the above equation, and is known from the relationship between the roots and the coefficients
Further, the S500 includes a region growing algorithm, which includes the following steps:
s501, selecting a seed point p (x)0，y0) The stack represents a seed area, and the seed point is pushed into the seed stack;
s502, stacking a first seed point pop in the seed stack, and traversing pixels in the neighborhood of the center 8 by taking the point as the center;
s503, judging whether the traversal pixel point is in the seed region, if not, judging whether the traversal pixel point meets the similarity of adjacent seed points, and if the pixel point (x, y) meets the similarity, pushing the pixel point (x, y) into a stack;
s504, repeating S502-S503 until the seed stack is empty.
Specifically, there are three elements that affect the region growing algorithm: selecting seed points; selecting a search path; and (5) judging the pixel similarity.
Selecting seed points: in general, the region growing algorithm is a semi-interactive segmentation algorithm, and requires a user to select a seed point, or a seed point calculated by other algorithms.
For the selection of the search path: the search path generally selects adjacent pixels, and takes a two-dimensional image as an example, the search path is generally an 8-neighborhood search or a 4-neighborhood search; taking a three-dimensional image as an example, a 26-neighborhood search or a 6-neighborhood search is common.
And (3) judging the pixel similarity: the similarity generally uses the degree of similarity of pixel values as a judgment standard, for example, a certain gray scale range can be set as a similarity standard; it is also possible to calculate that a certain shape or property is satisfied as the criterion.
Further, in step S500, for the triangle mesh hole filling method of Poisson 'S equation, the algorithm first needs to generate an initialized hole filling mesh according to the hole boundary, and then modifies the geometry of the triangle patch in the hole filling mesh through normal estimation and Poisson' S equation, so that it can adapt to and fuse with the surrounding original mesh. The algorithm mainly comprises the following steps: as shown in the figure, the hole boundary is detected and the hole filling mesh is initialized, and because the initialized hole filling mesh cannot be effectively fused with the mesh around the original hole, the vertex position of the hole filling mesh needs to be adjusted so that the hole filling mesh and the original mesh are in smooth transition; calculating an expected normal direction of vertexes in the hole-filling grid, taking the normal direction of an original grid hole boundary as the normal direction of the hole-filling grid boundary, constructing a Laplace equation delta f which is 0, and solving the normal distribution of vertexes inside the hole-filling grid, wherein if f represents a scalar on each vertex, a Laplace operator at the vertex xi on a grid domain is defined as follows (without considering the area influence):
where N1(xi) represents the 1-ring neighborhood point of vertex xi, αijAnd βijIs an edge eij2 corresponding diagonal corners; based on the expected normal direction, rotating the triangular patches in the hole filling grid, calculating to obtain the expected normal direction of the vertexes in the hole filling grid, and then further obtaining the expected normal direction of the triangular patches, wherein the expected normal direction of the triangular patches is the average value of the expected normal directions of the three vertexes, and then all the triangular patches in the hole filling grid rotate according to the expected normal direction, and the rotation parameter calculation method comprises the following steps: as shown, assume ni、ni’And ciFor a triangular patch fiOriginal normal, desired normal and center of gravity position, niAnd ni’The cross direction a of the triangle patch fiIn the direction of the axis of rotation of (c), niAnd ni’The included angle phi between the two is a triangular patch fiIs rotated, then the triangular patch fiWill be given by ciRotating the rotating shaft a to a new position by a rotation angle phi as a rotation center; adjusting the vertex position of the hole-filling grid based on a Poisson equation, tearing the hole-filling grid by rotating a triangular surface sheet of the hole-filling grid, reconstructing the hole-filling grid into a continuous grid curved surface by using the Poisson equation, calculating a gradient field of the torn grid when the Poisson equation is established, taking the gradient field as a guide field of the Poisson equation, and adjusting the vertex position of the grid,
f is the position of the peak of the adjusted grid to be solved, and w is the gradient field of the torn grid; gradient operators, assuming f denotes a scalar on each vertex, then the gradient operator of the scalar field f within any triangular patch T on the mesh domain is defined as follows:
wherein the expression for the gradient of the basis function ▽ Φ i is
± denotes the vector is rotated 90 degrees counter-clockwise, AT denotes the area of the triangle T, the divergence operator: assuming that w represents the vector at each triangle, the divergence operator of the vector field w at the vertex xi in the grid domain is defined as follows:
t1(xi) represents a 1-ring neighborhood triangle of the vertex xi, AT represents the area of the triangle T;
and S605, performing iterative region growing based on the curvature of the defect region obtained in the S400 to obtain a 3D point cloud patch. Specifically, as shown in fig. 1, iterative region growing is performed based on the optimal contour curvature obtained in step S400 to obtain an optimal 3D point cloud patch, where the point cloud patch generation formula is
Wherein, the area of A (v) depends on the size of the triangle vertex angle, and theta represents the triangle vertex angle of the calculated area.
Claims (10)
1. A skull-based CBCT image restoration guide plate generation method is characterized by comprising the following steps:
s100, obtaining a three-dimensional skull model through optical scanning, and establishing a projection image and a CBCT file sequence;
s200, analyzing a CBCT file sequence through a Marching Cube and a region growing algorithm;
s300, detecting a skull defect edge point cloud by an OpenCV edge extraction method, and generating a mirror image by using symmetrical points in the point cloud;
s400, calculating the curvature of the defect area in the mirror image, wherein the curvature comprises an average curvature, a Gaussian curvature and a main curvature;
s500, executing a region growing algorithm based on the curvature in the S400 to form a first-stage repair guide plate;
s600, optimizing the edge of the first-stage repairing guide plate based on a Poisson equation hole-filling smoothing method to form a final repairing guide plate.
2. The skull-based CBCT image restoration guide generation method according to claim 1, wherein the Marching Cube analysis process comprises the following steps:
s201, reading the original data into a specified array after preprocessing;
s202, extracting a unit body from the grid data body to form a current unit body, and simultaneously acquiring all information of the unit body, wherein the information comprises a vertex function value of a boundary and a point cloud coordinate position of the unit body;
s203, comparing the function value of the vertex in the current unit body with the given isosurface value C to obtain a state table of the unit body;
s204, finding out the unit body edge intersected with the isosurface according to the current state table index of the unit body, and calculating the position coordinates of each intersection point by adopting a linear interpolation method;
s205, solving the normal vector of the vertex in the current unit body by using a center difference method, and obtaining the normal vector of each vertex of the triangular patch by using a linear interpolation method;
and S206, drawing the isosurface image according to the coordinates of the vertexes of the triangular patches and the vertex normal vectors.
3. The guide plate repairing method based on the CBCT image of the skull as recited in claim 1 or 2, wherein the step of executing the region growing algorithm in S200 comprises controlling the threshold value at 770-25584 interval, and reconstructing a 3D model to view the main anatomical structure of CBCT.
4. The guide plate repairing method based on the skull CBCT image according to claim 1 or 2, wherein the OpenCV edge extraction method in S300 comprises the following steps:
generating a group of normalized Gaussian kernels through a discretized Gaussian function, and then carrying out weighted summation on each point of a gray matrix in the image data based on the Gaussian kernels;
highlighting points with significant changes in the intensity values of the neighborhood of the gray points in the image data through an enhancement algorithm;
and eliminating the noise edge points by a thresholding method.
5. The guide plate repairing method based on the CBCT image of the skull as recited in claim 1 or 2, wherein the mirror image generation of the symmetry points in S300 comprises the extraction of feature points, and the extraction formula is as follows:
wherein p isi,pjAre two points of symmetry, hi,hjIs the result after normalization, Ci,jIs the vector obtained, C is extractedi,jMaximum value to obtain feature points.
6. The guide plate repairing method based on the CBCT image of the skull as recited in claim 5, wherein the mirror image generation of the symmetry points in S300 comprises the extraction of feature points, and the extraction formula is as follows:
wherein p isi,pjAre two points of symmetry, hi,hjIs the result after normalization, Ci,jIs the vector obtained, C is extractedi,jMaximum value to obtain feature points.
7. The guide plate repairing method based on the CBCT image of the skull as recited in claim 1 or 6, wherein the mean curvature of the defect region in S400 is obtained by the following formula:
wherein, P,Q is the curvature mean value of the characteristic point, C is the form distance, and T is TPS and then the transformation matrix weight is estimated.
8. The guide plate repairing method based on the CBCT image of the skull as recited in claim 5, wherein the mean curvature of the defect region in S400 is obtained by the following formula:
wherein, P and Q are curvature mean values of the characteristic points, C is a form distance, and T is TPS estimated transformation matrix weight.
9. The guide plate repairing method based on the CBCT image of the skull as recited in claim 1 or 2, wherein the S500 comprises a region growing algorithm, comprising the following steps:
s501, selecting a seed point p (x)0，y0) The stack represents a seed area, and the seed point is pushed into the seed stack;
s502, stacking a first seed point pop in the seed stack, and traversing pixels in the neighborhood of the center 8 by taking the point as the center;
s503, judging whether the traversal pixel point is in the seed region, if not, judging whether the traversal pixel point meets the similarity of adjacent seed points, and if the pixel point (x, y) meets the similarity, pushing the pixel point (x, y) into a stack;
s504, repeating S502-S503 until the seed stack is empty.
10. The skull-based CBCT image guide plate repairing method according to claim 1, wherein the step S600 comprises generating an initialized hole filling grid according to the hole boundaries in the triangular grid, and then modifying the geometric shape of the triangular patch in the hole filling grid by normal estimation and Poisson' S equation to make it able to adapt and fuse with the surrounding original grid, the process comprising the following steps:
s601, detecting hole boundaries and initializing hole filling grids;
s602, calculating an expected normal direction of the vertexes in the hole filling grid, constructing a Laplace equation delta f which is 0, and solving normal distribution of the vertexes inside the hole filling grid, wherein a Laplace operator is as follows:
wherein N is1(xi) Representing a vertex xi1 ring neighborhood point of αijAnd βijIs an edge eij2 corresponding diagonal corners;
s603, rotating the triangular patches in the hole filling grid based on the expected normal direction;
s604, adjusting the vertex position of the hole filling grid based on a Poisson equation, and calculating the gradient field of the three-dimensional grid, wherein the formula is as follows:
f is the position of the peak of the adjusted grid to be solved, and w is the gradient field of the torn grid;
the gradient operator is:
wherein the expression of the basis function gradient ▽ Φ i is:
it represents the vector rotated 90 degrees counterclockwise, AT represents the area of the triangle T;
divergence operator:
wherein T1(xi) represents the 1-ring neighborhood triangle of the vertex xi, and AT represents the area of the triangle T;
and S605, performing iterative region growing based on the curvature of the defect region obtained in the S400 to obtain a 3D point cloud patch.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911148810.0A CN110942433B (en) | 2019-11-21 | 2019-11-21 | Repairing guide plate generation method based on skull CBCT image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911148810.0A CN110942433B (en) | 2019-11-21 | 2019-11-21 | Repairing guide plate generation method based on skull CBCT image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110942433A true CN110942433A (en) | 2020-03-31 |
| CN110942433B CN110942433B (en) | 2023-11-03 |
Family
ID=69907020
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911148810.0A Active CN110942433B (en) | 2019-11-21 | 2019-11-21 | Repairing guide plate generation method based on skull CBCT image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110942433B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111768353A (en) * | 2020-06-30 | 2020-10-13 | 北京华捷艾米科技有限公司 | Hole filling method and device for three-dimensional model |
| CN114708973A (en) * | 2022-06-06 | 2022-07-05 | 首都医科大学附属北京友谊医院 | Method for evaluating human health and related product |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6246784B1 (en) * | 1997-08-19 | 2001-06-12 | The United States Of America As Represented By The Department Of Health And Human Services | Method for segmenting medical images and detecting surface anomalies in anatomical structures |
| CN1523530A (en) * | 2003-09-10 | 2004-08-25 | 北京工业大学 | Method for preparing titanium alloy skull dummy |
| CN102302803A (en) * | 2011-09-15 | 2012-01-04 | 周强 | Human skull repairing scaffold and preparation method thereof |
| CN105046750A (en) * | 2015-08-24 | 2015-11-11 | 杭州美齐科技有限公司 | Method for automatically segmenting whole dental triangular mesh model |
| CN106780591A (en) * | 2016-11-21 | 2017-05-31 | 北京师范大学 | A kind of craniofacial shape analysis and Facial restoration method based on the dense corresponding points cloud in cranium face |
| CN107874831A (en) * | 2017-11-21 | 2018-04-06 | 四川大学 | A kind of cranium Maxillary region guide plate design method based on implicit function |
| CN110189352A (en) * | 2019-05-21 | 2019-08-30 | 重庆布瑞斯科技有限公司 | A kind of root of the tooth extracting method based on oral cavity CBCT image |
| CN110223308A (en) * | 2019-04-15 | 2019-09-10 | 东南大学 | A kind of stack position point cloud localization method increased based on edge detection and region |
-
2019
- 2019-11-21 CN CN201911148810.0A patent/CN110942433B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6246784B1 (en) * | 1997-08-19 | 2001-06-12 | The United States Of America As Represented By The Department Of Health And Human Services | Method for segmenting medical images and detecting surface anomalies in anatomical structures |
| CN1523530A (en) * | 2003-09-10 | 2004-08-25 | 北京工业大学 | Method for preparing titanium alloy skull dummy |
| CN102302803A (en) * | 2011-09-15 | 2012-01-04 | 周强 | Human skull repairing scaffold and preparation method thereof |
| CN105046750A (en) * | 2015-08-24 | 2015-11-11 | 杭州美齐科技有限公司 | Method for automatically segmenting whole dental triangular mesh model |
| CN106780591A (en) * | 2016-11-21 | 2017-05-31 | 北京师范大学 | A kind of craniofacial shape analysis and Facial restoration method based on the dense corresponding points cloud in cranium face |
| CN107874831A (en) * | 2017-11-21 | 2018-04-06 | 四川大学 | A kind of cranium Maxillary region guide plate design method based on implicit function |
| CN110223308A (en) * | 2019-04-15 | 2019-09-10 | 东南大学 | A kind of stack position point cloud localization method increased based on edge detection and region |
| CN110189352A (en) * | 2019-05-21 | 2019-08-30 | 重庆布瑞斯科技有限公司 | A kind of root of the tooth extracting method based on oral cavity CBCT image |
Non-Patent Citations (8)
| Title |
|---|
| EZZET AL-HUJAZI等: "Range Image Segmentation with Applications to Robot Bin-Picking Using Vacuum Gripper" * |
| JUN WANG等: "A Morse-Theory Based Method for Segmentation of Triangulated Freeform Surfaces" * |
| WEI ZHAO等: "A robust hole-filling algorithm for triangular mesh" * |
| 庞正雅等: "一种改进的点云数据三维重建算法" * |
| 张琪: "三维网格模型层次分割及骨架提取" * |
| 禇国民等: "钛金属板颅骨修复体数字化制造研究" * |
| 赵伟: "自由形状特征的重用与抑制" * |
| 钟昌康等: "基于曲率的三角网格模型分割算法" * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111768353A (en) * | 2020-06-30 | 2020-10-13 | 北京华捷艾米科技有限公司 | Hole filling method and device for three-dimensional model |
| CN111768353B (en) * | 2020-06-30 | 2023-11-03 | 北京华捷艾米科技有限公司 | Hole filling method and device for three-dimensional model |
| CN114708973A (en) * | 2022-06-06 | 2022-07-05 | 首都医科大学附属北京友谊医院 | Method for evaluating human health and related product |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110942433B (en) | 2023-11-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110189352B (en) | Tooth root extraction method based on oral cavity CBCT image | |
| JP7493464B2 (en) | Automated canonical pose determination for 3D objects and 3D object registration using deep learning | |
| US8571278B2 (en) | System and methods for multi-object multi-surface segmentation | |
| US8913060B2 (en) | Systems and methods for extracting a curve-skeleton from a volumetric image of a vessel | |
| US9367958B2 (en) | Method and apparatus for correction of errors in surfaces | |
| Xia et al. | Individual tooth segmentation from CT images scanned with contacts of maxillary and mandible teeth | |
| US20040228529A1 (en) | Systems and methods for providing automatic 3D lesion segmentation and measurements | |
| Jones | Facial Reconstruction Using Volumetric Data. | |
| CN111968146B (en) | Three-dimensional dental mesh model segmentation method | |
| WO2020263987A1 (en) | Processing digital dental impression | |
| Chartrand et al. | Semi-automated liver CT segmentation using Laplacian meshes | |
| Chen et al. | Shape registration with learned deformations for 3D shape reconstruction from sparse and incomplete point clouds | |
| CN110942433B (en) | Repairing guide plate generation method based on skull CBCT image | |
| Kretschmer et al. | ADR-anatomy-driven reformation | |
| Jung et al. | Combining volumetric dental CT and optical scan data for teeth modeling | |
| CN111127488A (en) | Method for automatically constructing patient anatomical structure model based on statistical shape model | |
| Yao | A statistical bone density atlas and deformable medical image registration | |
| Liu et al. | Tracking-based deep learning method for temporomandibular joint segmentation | |
| Schoenemann et al. | Validation of plaster endocast morphology through 3D CT image analysis | |
| Zhao et al. | A frame of 3D printing data generation method extracted from CT data | |
| CN110689080A (en) | Planar atlas construction method of blood vessel structure image | |
| Abdolali et al. | Mandibular canal segmentation using 3D Active Appearance Models and shape context registration | |
| JP5954846B2 (en) | Shape data generation program, shape data generation method, and shape data generation apparatus | |
| Wei et al. | Automatic mesh fusion for dental crowns and roots in a computer-aided orthodontics system | |
| Song et al. | Confidence Surface-Based Fine Matching Between Dental Cone-Beam Computed Tomography Scan and Optical Surface Scan Data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |