A three-dimensional wound surface finite element development method for mapping skin material properties

Abstract

The present invention discloses a three-dimensional wound surface finite element development method for mapping skin material properties, which relates to the field of reverse engineering. The method includes four parts: topological relationship reconstruction of wound surface mesh model, hierarchical division of surface structure, isomorphic development of wound surface and optimization of unit material model. The method can realize the reading and processing of wound surface data, rapid topological reconstruction of the model, hierarchical division of wound surface structure, morphological development of wound surface model, iterative optimization based on unit material model and multi-resolution smoothing optimization of flattened contour, and realize the method of determining flap shape from wound surface morphology. The three-dimensional wound surface finite element development method of the present invention has universality and effectiveness for disordered and irregular three-dimensional wound surfaces, and the method is stable and the calculation process is simple. The clinical application results show that it can effectively reduce the preparation time of preoperative plans, improve the work efficiency of doctors, and improve the accuracy of preoperative plans and the reliability of evaluation.

Description

Technical Field

[0001] The invention relates to the field of reverse engineering, in particular to a three-dimensional wound surface finite element development method for mapping skin material properties. Background Art

[0002] Surface development is the process of mapping a spatial surface to a specific plane. In industrial design, engineers often determine the shape of the finished blank by surface development, which is convenient for subsequent processing and manufacturing. According to the complexity of the morphology, surfaces can be divided into simple surfaces and complex surfaces. Among them, most simple surfaces can be developed, that is, the Gaussian curvature on the surface is zero everywhere, and they can be equidistantly transformed to a certain plane during the development process, and the development result is accurate and unique. Complex surfaces will inevitably produce deformations such as area and angle when mapped to a certain plane. For example, the clinical wound case involved in the present invention is a typical complex non-developable surface. In order to achieve the approximate development of complex non-developable surfaces, in the existing surface development methods, according to their characteristics, they can be divided into two categories: geometric methods and mechanical methods. The geometric method is based on the geometric characteristics of the surface, such as side length, area, angle, etc. It is developed to a certain plane through given rules, such as ensuring that the front and rear side lengths of the grid remain unchanged. The mechanical method regards the deformation of the surface during the development process as the strain generated by the surface material, and calculates its shape after development by analyzing the stress and strain of the surface.

[0003] Regarding the study of geometric methods, since Manning used the isometric tree theory to analyze the problem of upper expansion and cutting in 1980, domestic and foreign scholars have proposed more surface expansion methods and applied them to different fields. Bennis realized the geometric expansion of the surface based on geodesic curvature. First, the surface was divided into equal parameter grids, and then the geodesic curvature of the surface was calculated and expanded. This method was applied in computer-aided shoe design, but because it did not use an adaptive partitioning algorithm, it was inefficient and prone to cumulative errors. In 1993, Parida proposed a geometric expansion method for composite laminates. First, the surface was divided into multiple triangular grids, and then the grids were expanded one by one to the plane. The direction of cracks and overlaps on the plane was adjusted to meet the expansion constraints, but more cracks were generated during the expansion process, and the cumulative errors were more obvious. In 2015, in order to improve the design method of fiber composite materials, Kotaro divided the three-dimensional sheet model into two surfaces and flattened it through the parameterization method of the cross vector field. This method can handle large deformations of thick plate materials in the fiber direction and thickness direction, but whether it conforms to the actual deformation of the actual material remains to be verified.

[0004] Regarding the study of mechanical methods, in 1991, Shinmada used dynamic programming to approximate the transformation of a space surface to a plane. In the process, the surface to be developed was first divided into multiple strip regions based on the finite element method, and then each region was unfolded in turn. This method is not suitable for processing more complex surfaces. In 2001, Zhang Qilin derived an incremental finite element method for surface unfolding, which transformed the unfolding problem into the problem of finding a plane triangular unit set with a size closest to the space grid, and established the closest judgment criterion, that is, the residual node force in the triangular unit is zero, but the efficiency of this method in processing a large number of mesh nodes still needs to be improved.

[0005] In the current clinical preoperative flap design, the traditional sample cloth tracing method is relatively rough and takes a long time to prepare. The plan formulation process is more dependent on the doctor's clinical experience and habits, and the matching degree between the cut flap and the wound surface is poor. Therefore, considering introducing digital technology into the preoperative flap plan design, a three-dimensional wound surface finite element development method that maps the skin material properties is proposed. Summary of the invention

[0006] In view of the above problems and technical requirements, the inventors have proposed a three-dimensional wound surface finite element development method for mapping skin material properties. The technical solution of the present invention is as follows:

[0007] A three-dimensional wound surface finite element development method for mapping skin material properties comprises the following steps:

[0008] (1) Reading and preprocessing of wound surface mesh model data

[0009] In the clinic, laser scanners are used to collect wound surface data, and combined with the doctor's guidance, the wound surface mesh model required for preoperative flap design is extracted to generate the corresponding "*.stl" file. The point-edge-face class structure is designed to build a storage framework for surface information, and the topological association between points, edges and faces on the wound surface model is reconstructed using the hash storage structure. The hierarchical division of the surface structure is realized based on the hierarchical value, and the generalized hierarchical value is applied to cluster the mesh surface for subsequent iterative optimization of the surface vertex position.

[0010] (2) Morphological development of wound surface mesh model

[0011] For different types of wounds, the principle of selecting the base point for their development is formulated according to the level value. Based on the equal deformation method, the deformation of the material after the wound surface is developed is regarded as uniform and continuous, and the deformation caused by the interaction between the units on the surface during the development process is controlled.

[0012] Due to the different causes of wounds, such as external force, electrical injury, chemical burns, etc., their surface morphology is special. In the wound mesh model composed of triangular meshes, the equal deformation method has regional mesh deformation caused by the position mutation of individual vertices, which produces a larger cumulative error. Such vertices are defined as morphological distortion points. The unfolding of morphological distortion points takes the overall optimal unfolding of the wound as the starting point, uniformly deforms the adjacent neighborhood angles to amortize the mutations in the unfolding process, and based on the principle of adjacent edge invariance, ensures that the length of the neighborhood edges around the morphological distortion points remains unchanged. In order to reasonably and scientifically analyze and evaluate the unfolding results of the wound, from the perspective of geometric characteristics, the mean area and edge length errors between the meshes before and after unfolding are used to evaluate the flattening results of the wound.

[0013] (3) Optimization of wound surface development results based on unit material model

[0014] The unit material model is used to optimize the initial results of wound surface expansion. The three nodes of the triangular mesh unit are used to define the deformation displacement of each unit before and after expansion. The unit material model related to the deformation displacement of the unit node on the plane mesh is established based on the material matrix. The wound flattening shape optimization model with the total deformation energy of the model as the objective function is established according to the equivalent node load. The Gaussian elimination method is used to solve the optimized displacement of each unit node on the wound surface. The two-dimensional equal bandwidth storage method is used to improve the solution accuracy and calculation efficiency of the wound surface optimization model during the solution process. The contour curve of the preoperative flap design is fitted based on the cubic non-uniform rational B-spline (NURBS) curve interpolation method to solve the problem of uneven distribution of the skin contour point set after expansion, and obtain a smooth and fair boundary curve suitable for medical scenarios.

[0015] The beneficial technical effects of the present invention are:

[0016] By establishing the topological association between points, edges and surfaces on the wound mesh model, and then dividing the mesh levels to enrich the expression of the model structure, the wound model morphology is unfolded to a given plane based on the equal deformation method, and the unit material model and multi-resolution smoothing method are used to achieve boundary fitting and optimization of the wound unfolding result. The final surface unfolding result can better match the skin wound. The three-dimensional wound finite element unfolding method of the present invention is universal and effective for disordered and irregular three-dimensional wounds, and the method is stable and the calculation process is simple. The clinical application results show that it can effectively reduce the preparation time of preoperative plans, improve the work efficiency of doctors, and improve the accuracy of preoperative plans and the reliability of evaluation. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a flow chart of the three-dimensional wound surface finite element development method provided in this application.

[0018] Figure 2 It is a schematic diagram of a triangular face group provided in this application.

[0019] Figure 3 This is a schematic diagram of the point-edge-face topology information storage method when reading a model file provided by this application.

[0020] Figure 4 This is a schematic diagram after the triangular mesh level division provided by this application is completed.

[0021] Figure 5 (a) is the spatial grid before expansion provided by the present application.

[0022] Figure 5 (b) is the unfolded plane grid provided by this application.

[0023] Figure 6 It is a schematic diagram of the three-node displacement of the triangular mesh unit provided in this application.

[0024] Figure 7 It is the unfolding result of the wrist wound surface morphology obtained by the equal deformation method provided by the present application.

[0025] Figure 8 This is the result of wrist wound expansion optimized by the unit material model provided in this application.

[0026] Fig. 9 It is the skin contour after boundary extraction provided by this application.

[0027] Fig.10 This is the skin contour provided by the present application after boundary fitting and multi-resolution smoothing.

[0028] Fig.11 It is a schematic diagram of the effect of skin contour fitting wound surface provided by the present application. DETAILED DESCRIPTION

[0029] The specific implementation of the present invention will be further described below in conjunction with the accompanying drawings.

[0030] This application discloses a three-dimensional wound surface finite element development method for mapping skin material properties, such as Figure 1 As shown, the specific steps include:

[0031] Step 1: Establish the topological relationship between points, edges and surfaces of the wound surface mesh model based on the wound surface data.

[0032] The wound surface model data extracted from the patient's wound is stored in the "*.stl" file, but the file can only record the direction and position information of each facet in an unordered manner, and the topological information between the faces on the original model has been lost. It is necessary to use the hash storage structure to reconstruct the topological association between points, edges and faces on the wound surface mesh model. This application takes the triangular facet group that is commonly found in the mesh model as an example, such as Figure 2The reconstruction process of topological association is implemented by compiling the function "Topo.cpp" in Visual Studio, as described below:

[0033] Define the point table CPoint, the patch vertex table CVertex, the edge table CEdge, the face table CTriangle and the topological hash function HashNode.

[0034] 1) Points are the basic elements that constitute the surface. In the process of reading STL files, the geometric information of the wound surface data containing points is stored in the CPoint class, and the corresponding calculation methods are added.

[0035] 2) The patch vertex table CVertex is required to not only retrieve the coordinate information of the corresponding point, but also record the topological relationship between the point and the surrounding points-edges-faces. Therefore, the adjacent point AdjPoint, the neighboring edge NeiEdge, the adjacent edge AdjEdge and the neighboring triangle NeiTri are introduced into the patch vertex table to enrich the topological association of the vertex.

[0036] 3) Edges are elements abstracted from the "*.stl" file. In order to enrich the topological information of the surface and meet the requirements of calculation methods such as area and curvature, the two vertex indices p[2], the neighborhood point NeiPoint and the adjacent triangle AdjTri index that make up the edge are added to the edge table CEdge.

[0037] 4) The "*.stl" file uses the face as the storage unit, and records the normal vector coordinates of each face, as well as the topological information such as the points and edges that make up the face. Therefore, the face normal vector Tnormal, three vertex indices p[3], three edge indices e[3], and adjacent triangle AdjTri indices are added to the face table CTriangle.

[0038] 5) The establishment of topological hash function HashNode uses hash function f v Convert the key value K of the point into the corresponding hash table address (i.e., hash table subscript) to access the point index in the patch vertex table. The conversion formula is:

[0039]

[0040] f v =KmodM

[0041] Among them, x, y, z are the coordinate components of the vertex in each direction, and M is the length of the hash table, which is a prime number. Similarly, using the hash function f eCalculate the key value and hash address of the edge and store them in the edge table CEdge. Read the triangle patch information in the "*.stl" file in sequence, store the normal vector and vertex coordinate information in the surface table CTriangle, and realize the rapid topological reconstruction of the wound surface mesh model.

[0042] The specific storage format is shown in Table 1. When describing the structural relationship of the wound mesh model, only the vertex coordinate data is stored and placed in the point table CPoint. The corresponding position information of the edge and triangle is constructed using the vertex index number to reduce computer storage space. The structure relationship of point-edge-face is as follows: Figure 3 shown.

[0043] Table 1 Recording method of topological relationship

[0044]

[0045] Step 2: Divide the wound surface mesh model into mesh levels based on the level value and determine the initial expansion base point.

[0046] The concept of level value comes from the description of the progressive change of the level of the surface, starting from the outer boundary of the surface and gradually approaching inward. Each point is assigned a level value vl, and points with the same level value are located at the same level. The corresponding level contour on the mesh surface is extracted according to the level value. Among them, the edges of two adjacent triangles are internal edges, and the edges of only one adjacent triangle are boundary edges. The two adjacent points of the boundary edge are boundary points, and the boundary triangle containing the boundary edge is a boundary triangle. The level value is used as the basis to realize the assignment and division of the grid points.

[0047] After traversing all points on the wound mesh surface, the maximum level value is obtained as vl max , re-assign the level values ​​of all vertices, the assignment formula is: vl = vl max -vl+1(2)

[0048] Define the grid points with the same level value vl as the same level. Therefore, the level contour in the grid surface can be extracted as follows Figure 4 As shown. The first layer is the layer contour composed of all points with a layer value of 1. Similarly, the points with a layer value of 4 are the boundary points of the surface, and the polylines formed by the boundary points are the boundary contours of the surface.

[0049] When mapping a spatial surface to a given plane, the initial expansion base point is the first point mapped to the plane, and the point with vl = 1 is taken as the initial expansion base point of the wound mesh model. When this point is not unique, the center position of the entire surface mesh is calculated, and the point closest to the center position is selected as the initial expansion base point.

[0050] Step 3: Set the clustering threshold ε to determine whether the surface needs further clustering.

[0051] In differential geometry, the method for determining a developable surface is that the Gaussian curvature of each point on it is always zero, but this is not enough to evaluate the developability of the surface. Generally speaking, the larger the discrete Gaussian curvature of the surface vertex, the more likely it is that the mesh will deform greatly after being unfolded to a plane. This application uses the following formula to describe the developability of the surface:

[0052]

[0053] Where A is the total area of ​​the mesh surface, A sum (P i ) represents point P i The total area of ​​all neighboring triangles, θ sum (P i ) is point P i The sum of the interior angles formed by all its neighboring triangles. i When it is a boundary point, take θ sum (P i )=2π. D(F) is between 0 and 1. When D(F) is closer to 0, it means that the surface is less developable and more likely to deform when unfolded.

[0054] The clustering threshold ε is set. When D(F) < ε, the surface is considered to have poor developability and no clustering is required. When D(F) ≥ ε, the surface is considered to have good developability and further clustering is required. The clustering threshold is based on the vertex Gaussian curvature value K in the wound surface mesh model. p The value range of is given by:

[0055]

[0056] Where θ is the sum of the interior angles of point P.

[0057] The clustering method uses K-means clustering, randomly selecting K points in the grid as the initial center of each class, integrating the remaining vertices into the nearest class according to their similarity with each cluster center, and iteratively updating the value of each cluster center until the expected clustering result is met.

[0058] Step 4: Starting from the initial expansion base point, based on the equal deformation method, the adjacent grid points are expanded layer by layer from the inside to the outside to the given plane to obtain the flattened shape of the wound surface.

[0059] The concept of iso-deformation originates from the assumption of uniform continuity of material deformation. In order to reproduce the actual situation of real materials in deformation such as stretching and simple bending, when the surface is divided into several triangular units by mesh, the deformation of the units when unfolding from space to plane is still uniform and continuous. The essence of the iso-deformation method is still the geometric method in the theory of surface unfolding, but its unfolding process conforms to the characteristics of uniform continuity of local areas of the surface in actual deformation. For any three-dimensional surface, a plane with topological equivalent can be found in two-dimensional space.

[0060] like Figure 5 As shown in (a), ΔA 0 B 0 C 0 With ΔB 0 C 0 D 0 Point A is two adjacent triangular mesh units in the mesh before expansion. 0 With point D 0 For edge B 0 C 0 Two neighboring points of ; Figure 5 (b) is the corresponding plane mesh after being unfolded by the equal deformation method, where ΔABC and ΔBCD are two topologically equivalent triangular mesh units.

[0061] On the spatial grid, record ΔB 0 C 0 D 0 , ΔA 0 D 0 C 0 , ΔA 0 B 0 D 0 , ΔA 0 B 0 C 0 The area is S A0 , S B0 , S C0 , S D0 , which can be calculated by the mesh edge length. The corresponding areas of ΔBCD, ΔADC, ΔABD, and ΔABC after equal deformation expansion are S A , S B , S C , S D , represented by the coordinates of point D. Based on the assumption of equal deformation, the unit area change rate before and after expansion is defined as τ, then:

[0062]

[0063] Among them, S i0 is the area of ​​the triangle formed by the three vertices except i before unfolding, S iis the area of ​​the unfolded triangle formed by the three vertices except i; the larger τ is, the greater the deformation of the corresponding unit after the mesh is unfolded.

[0064] After the expansion, the area of ​​the triangular mesh unit can also be expressed by the corresponding three vertex coordinates (x i ,y i )express:

[0065]

[0066] Among them, vertex D is the adjacent point to be expanded on the plane, and the coordinates (x D ,y D )Unknown; T A , T B , T C They represent the areas of ΔBCD, ΔADC, and ΔABD after expansion, respectively.

[0067] Formula (5) and Formula (6) are both used as methods to solve the unit area of ​​the triangular mesh, so the unit areas obtained are equal, which leads to:

[0068] S i =T i i=A,B,C (7)

[0069] The equation of formula (7) belongs to an overdetermined system of equations. In order to calculate its optimal approximate solution, the distance metric M is introduced:

[0070] M(S,T)=(S A -T A ) 2 +(S B -T B ) 2 +(S C -T C ) 2 (8)

[0071] Through M, the unit area obtained by equation (7) and equation (8) is closest, that is, the minimum distance between S and T under the distance metric M is solved. It can be seen from equation (8) that:

[0072] Since M is non-negative, we solve the extreme point of M to determine the minimum value of M, and combine it with formula (5) to get the unknown number x D ,y D The expression is:

[0073]

[0074] Among them, vertex D 0 are the adjacent points to be developed on the surface; so far, the points to be developed on the surface D are determined one by one through the equal deformation method 0The plane coordinates (x D ,y D ) to realize the planar unfolding of the mesh surface.

[0075] Optionally, the above step 4 further includes processing of morphological distortion points, specifically:

[0076] The concept of developability is introduced by combining the uniform continuity assumption of the equal deformation method. The developability function of the vertices in the triangular mesh unit is:

[0077]

[0078] Among them, θ i It refers to the inner angle formed by the i-th neighborhood triangle of a vertex and the vertex. A threshold η is set. When D(v)≤η, it means that the point has not undergone a mutation and is a regular point. When D(v)>η, it means that the point has undergone a mutation and is defined as a morphological distortion point, but does not include boundary points.

[0079] When processing the expansion of morphological deformity points on the wound surface, the overall optimal expansion of the wound surface is taken as the starting point, and the expansion mutation of the morphological deformity points is evenly spread to its adjacent points in the form of neighborhood angles. The idea of ​​equal deformation is used to effectively reduce the degree of deformation in the local area, but when the expandability of the morphological deformity point is worse, the deformation reduced by spreading is still large. Therefore, the present application considers adding the adjacent edge invariance principle on the basis of the equal deformation method, and placing the morphological deformity point at the center of the neighborhood polygon formed by its neighboring points while ensuring that the length of each neighborhood edge of the morphological deformity point remains unchanged during expansion, and keeping the angle ratio formed by each adjacent edge of the morphological deformity point unchanged.

[0080] Optionally, the above step 4 further includes a reasonable and accurate analysis and evaluation of the expansion results, specifically:

[0081] The use of equal deformation method can avoid the appearance of flipped triangles as much as possible. The hierarchical structure established based on the hierarchical value allows the deformation of the mesh to be evenly distributed after unfolding, and obtains an unfolding result with good visual effect. However, a scientific and reasonable judgment standard is still needed to evaluate the quality of the surface unfolding results. This application establishes an evaluation standard for surface unfolding based on the change in the area or side length of the triangular mesh unit to achieve error analysis of the unfolding results.

[0082] The surface unfolding error is the mean of the errors corresponding to each vertex on the surface. The edge length error of the vertex refers to the length change error of all adjacent edges of the point before and after unfolding. The mean edge length error of the surface after unfolding is:

[0083]

[0084] Where m′ is the number of vertices on the surface mesh, n is the number of adjacent edges of point i, and L jis the length of the jth edge corresponding to point i in the unfolded plane mesh, L j 0 is the length of the j-th edge corresponding to point i in the spatial grid before expansion.

[0085] The area error of a vertex refers to the area change error of all the neighborhood triangles of the point before and after expansion. The mean area error of the surface after expansion is:

[0086]

[0087] Among them, k is the number of triangles in the neighborhood of point i, S j and S j 0 are the areas of the j-th neighborhood triangle corresponding to point i in the plane mesh after expansion and the space mesh before expansion, respectively.

[0088] Step 5: Establish a unit material model based on the unit node deformation and displacement of the wound surface mesh model.

[0089] The approximate unfolding of developable surfaces is inevitably accompanied by the deformation of mesh units. What standard should be used to measure unit deformation is the first issue that needs to be considered before optimizing the unfolding results. Traditional methods usually describe the overall mesh deformation based on the change in the length of the mesh unit side, such as the spring mass model. This measurement method simply adds the force model to the geometric deformation and cannot more realistically represent the deformation of the surface during actual unfolding. In the finite element method, since the three nodes of the triangular mesh unit have a strong adaptability to the boundary of the elastic body, the stress distribution law inside the elastic body is often reflected through the triangular unit when applied, such as Figure 6 As shown in the figure, the nodes of the triangular mesh unit are numbered counterclockwise (in the order of A, B, C). Each triangular mesh unit in the model has 6 node displacement components. The node displacement array of each triangular mesh unit after expansion is expressed as:

[0090]

[0091] Among them, δ i (i=A, B, C) represents the displacement of three nodes in the triangular mesh unit, which is composed of the displacement components (u i ,v i )express.

[0092] According to the displacement of each node in the triangular mesh unit and the corresponding node coordinates, the relationship between the displacement of any point in the unit and the node displacement is determined. The matrix form of the unit deformation displacement after expansion is:

[0093]

[0094] Among them, E is the second-order identity matrix, N A 、NB 、N C is a 2×2 shape function matrix, N is a 2×6 shape function matrix, then the expansion strain of the triangular mesh element is:

[0095]

[0096] Among them, ε x , ε y , γ xy is the normal strain and shear strain of the unit in the x and y directions; B is the strain matrix of the triangular mesh unit, and B = LN, L is the differential operator in the plane problem. B is only related to the node coordinates after its expansion. As long as the plane coordinates of the mesh unit on the model are determined after expansion, the strain matrix of the current unit can be obtained by establishing the relationship between the unit node displacement and its plane coordinates. This matrix is ​​also called a constant matrix. It is difficult to accurately reflect its changes in areas with large stress gradients. The number of units and nodes can be appropriately increased.

[0097] In the plane problem, the total potential energy functional established according to the minimum potential energy principle is expressed as:

[0098]

[0099] Among them, Q is the total potential energy of the system, which is the sum of the deformation potential energy of the elastic body and the potential energy of the external force. In the three-dimensional wound expansion problem, the model needs to be expanded to a given plane without external force, so the volume force f and the area force T on the boundary are both zero. The total potential energy of the original elastic body only includes the deformation energy, so the total potential energy functional of the triangular mesh unit, that is, the unit deformation energy, is:

[0100]

[0101] Among them, K e is the stiffness matrix of the triangular mesh unit; D is the elastic coefficient matrix of the triangular mesh unit. Due to the complex material properties of soft tissues such as skin, it is difficult to express them with a simple matrix. Here, it is first simplified to an ideal elastic model. This is reasonable because the process of optimizing the initial unfolding result is not the actual deformation process of the material during unfolding, but only looking for an optimal wound unfolding state, without paying attention to the deformation history of the material. Therefore, in the process of unfolding from space to plane, the mechanical performance of the triangular unit on the surface is considered to be isotropic, and the basic expression of its elastic coefficient matrix is:

[0102]

[0103] As a transversely anisotropic material, skin has three mutually perpendicular symmetry planes. The fiber arrangement direction and the direction perpendicular to the fiber distribution are defined as the main directions of the material. The matrix form of its generalized Hooke's law is:

[0104]

[0105] Since the wound surface optimization process is simplified to the solution process of a two-dimensional plane problem, combined with formula (19), the elastic coefficient matrix of the skin material under transverse anisotropy conditions is obtained as follows:

[0106]

[0107] Among them, G xy is the shear modulus between the two material directions x and y; E x 、E y Respectively refers to the elastic modulus of the material along the x-axis and y-axis; μ xy Refers to the Poisson's ratio in the y-axis direction when tension is applied in the x-axis direction, μ yx Similarly, the elastic modulus and Poisson's ratio satisfy the following relationship:

[0108] E x μ yx =E y μ xy (twenty two)

[0109] In fact, when dealing with the problem of wound surface optimization, there is a situation where the main direction of the material is inconsistent with the direction of the coordinate axis. At this time, it is necessary to convert the stress on the coordinate axis and the material axis. The expression is:

[0110]

[0111] Among them, σ x , σ y are the normal stresses of the unit in the x and y directions, τ xy is the shear stress on the element in the xy plane.

[0112]

[0113] Similarly, the strain conversion relationship between the main direction of the material and the coordinate axis direction is:

[0114]

[0115] Substituting formula (25) back into formula (23), we obtain:

[0116]

[0117] in,

[0118]

[0119] Remember D 0 =T 1 (θ)DA -1 T1 Substituting (-θ)A into equation (18), the deformation energy of the triangular mesh unit is obtained as follows:

[0120]

[0121] The total deformation energy of the wound surface mesh model is the sum of the deformation energies of each triangular mesh unit, that is:

[0122]

[0123] Based on formula (28), the deformation potential energy distribution model of the unit and the whole after the three-dimensional wound surface is unfolded is established, which is the basis for realizing the optimization of the wound surface unfolding shape with deformation energy as the objective function.

[0124] Step 6: According to the equivalent node load, a wound flattening shape optimization model is established with the total deformation energy of the unit material model as the objective function, and the optimized displacement of each unit node on the wound surface is solved.

[0125] During the unfolding process of the wound surface model, the stress and strain calculated based on the deformation and displacement of the mesh unit nodes are regarded as a load acting between the nodes on the unfolding plane, that is, the unit equivalent node load. In order to describe the loading of the equivalent node load on the mesh unit nodes before and after unfolding, the general form of the potential energy functional of the triangular mesh unit is established as:

[0126]

[0127] Among them, P e is the unit equivalent node load matrix, and P e =K e δ e In fact, after the wound mesh model is unfolded, the unit deformation potential energy obtained by the equal deformation method is not necessarily the minimum value, that is, it is not necessarily the best unfolding result. If the node deformation displacement of each unit after unfolding can be properly adjusted to reduce the overall deformation potential energy and tend to the minimum stable value, it means that the unfolding shape is optimized in a good direction, which can improve the unfolding quality of the wound model and improve the design method of the flap shape.

[0128] Therefore, a new unit node displacement is introduced into the wound surface mesh model:

[0129] δ=[u 1 v 1 … u i v i … u n v n ] T (30)

[0130] Then the overall deformation energy of the model after expansion is:

[0131]

[0132] Among them, G is the transformation matrix, which realizes the unit stiffness matrix K e The same order integration with the global stiffness matrix K is expressed as:

[0133]

[0134] According to the variational principle, the condition for the functional to take an extreme value is that its first variation is 0, so we get:

[0135] Kδ=λP (32)

[0136] in,

[0137]

[0138] λ is the adjustment amplitude of displacement, which is used to adjust the number of model optimizations. A small adjustment amplitude means that the node displacement is small, and the number of optimizations increases. A large adjustment amplitude means that the node displacement is large, and the number of optimizations decreases. However, too large an adjustment amplitude may cause distortion of the final optimization result. Therefore, it is necessary to reasonably select the λ value, such as 0<λ<0.5.

[0139] Finally, a suitable λ is substituted into formula (32) to solve the optimized displacement δ of each grid unit node on the wound surface.

[0140] Step 7: Fit the optimized displacement of each unit node on the wound surface to obtain the optimized wound flattening boundary.

[0141] The contour points on the optimized wound surface unfolding shape are the type value points of the unfolded skin contour curve, which are back-calculated through the non-uniform rational B-spline control points; the first and last type value points are determined by the method of heavy nodes, and the remaining contour points are used as segmentation points in the fitting curve. The shape of the curve between adjacent type value points is described in the form of parametric equations; the skin contour curve formed by the optimized displacement of each grid unit node is defined as the original curve S m , where m is the scale and the node vector is denoted as U = [0,0,0,0,u 4 ,…,u m+2 ,1,1,1,1], the adjacent type value points are regarded as the minimum segment in the B-spline curve, and the distance between the adjacent type value points determines the distribution of the node vector. The standard accumulated chord length method is used to determine the next node vector. The calculation method is:

[0142]

[0143] Among them, q i Represented as data in a node vector.

[0144] According to the NURBS basis function of each segment of the curve, the relationship between the control point and the shape value point is:

[0145]

[0146] Among them, u i+3 is the starting node, i∈[1,m-1]; C i (i=0,1,2,…,m+2) is the corresponding i-th control point, V i is the data point, B represents the spline basis function; since the skin wound contour is a typical closed curve, the first and last points are the same, that is, V 0 =V m , a total of m-1 equations; the first and last three control points overlap in sequence, and there is C m =C 0 , C m+1 =C 1 , C m+2 =C 2 , the number of unknowns also becomes m-1; the number of equations is equal to the number of unknowns, so all the control points of the skin contour curve are obtained by solving the following equation, and the specific shape fitting of the curve is realized by using the obtained control points to obtain the optimized wound flattening boundary;

[0147]

[0148] Where P (i＝0,...,m-1) It is the type value point.

[0149] Step 8: Use a multi-resolution smoothing algorithm to smooth the wound surface boundary.

[0150] Construct a low-pass filter group to filter out the high-frequency part of the original curve to obtain a smooth low-frequency curve. The algorithm is: n =[<Φ n |Φ n >] -1 [<Φ n |Φ m >]C m ;

[0151] Among them, [<Φ n |Φ n >] -1 [<Φ n |Φ m >] A low-pass filter group constructed to achieve curve smoothing, C m , C n are the control points of the target curve before and after smoothing respectively; the control points are substituted into equation (35) to complete the drawing of the smoothed curve, namely:

[0152]

[0153] Step 9: Use the smooth wound surface to flatten the border and formulate the preoperative flap design.

[0154] The effectiveness of the method proposed by the present invention is illustrated below in combination with specific test data:

[0155] contrast Figure 7 and Figure 8 It can be seen that the comparison of the data before and after optimization of the wound surface expansion result is shown in Table 2. The experimental data proves that the error analysis of the optimization result of the unit material model obtained by the method provided in this application is smaller, and the above method is feasible. In most cases, the boundary contours of the morphological expansion results of the three-dimensional wound surface have shortcomings such as not being smooth enough and not easy to perform intraoperative cutting. Affected by such factors, it is necessary to further optimize the contour of the wound surface expansion result before the actual flap is cut, in order to fit a contour curve that is more suitable for preoperative flap design in medical scenarios. Fig. 9 The skin contour of the wrist wound surface after boundary extraction is shown in Figure 1. The skin contour is obtained by boundary fitting and multi-resolution smoothing. Fig.10 As shown, the schematic diagram of the contour and the corresponding skin wound fitting effect is as follows Fig.11 As shown in the figure, the surface unfolding result after surface optimization can better fit the skin wound surface.

[0156] Table 2 Wrist wound development results data

[0157]

[0158] In summary, by establishing the topological association between points, edges and surfaces on the wound mesh model, the mesh level is divided to enrich the expression of the model structure, the wound model morphology is unfolded to a given plane based on the equal deformation method, and the unit material model and multi-resolution smoothing method are used to achieve boundary fitting and optimization of the wound unfolding result. The final surface unfolding result can better match the skin wound. The three-dimensional wound finite element unfolding method of the present invention is universal and effective for disordered and irregular three-dimensional wounds, and the method is stable and the calculation process is simple. The clinical application results show that it can effectively reduce the preparation time of the preoperative plan, improve the work efficiency of doctors, and improve the accuracy of the preoperative plan and the reliability of the evaluation.

[0159] The above is only a preferred embodiment of the present application, and the present invention is not limited to the above embodiments. It is understood that other improvements and changes directly derived or associated by those skilled in the art without departing from the spirit and concept of the present invention should be considered to be included in the protection scope of the present invention.

Claims

1. A three-dimensional wound surface finite element development method for mapping skin material properties, characterized in that: The method comprises: Establishing a topological relationship between points, edges and faces of a wound surface mesh model based on wound surface data; wherein the wound surface data includes direction and position information of each face piece; Divide the wound surface mesh model into mesh hierarchies based on the hierarchical value to determine an initial expansion base point; Starting from the initial expansion base point, adjacent grid points are expanded layer by layer from the inside to the outside to a given plane based on the equal deformation method to obtain a flattened shape of the wound surface; A unit material model is established based on the deformation displacement of the unit nodes of the wound mesh model, and a wound flattening shape optimization model is established based on the equivalent node load, with the total deformation energy of the unit material model as the objective function, to solve the optimized displacement of each unit node on the wound surface, and the optimized wound flattening boundary is obtained through fitting; wherein the equivalent node load is a load acting between each node on the unfolded plane, which is obtained by calculating the stress and strain calculated based on the deformation displacement of the unit node; The method of establishing a unit material model based on the unit node deformation displacement of the wound surface mesh model includes: The mesh on the wound surface mesh model is regarded as a triangular mesh unit, and the node displacement array of each triangular mesh unit after expansion is expressed as: Among them, δ i (i=A, B, C) represents the displacement of the three nodes in the triangular mesh unit, which is represented by the displacement component (u i ,v i )express; When the wound surface is unfolded to a plane, assuming that the model is unfolded to a given plane without external force, the total potential energy of the original elastic body system only includes the deformation energy, and the deformation energy of the triangular mesh unit is: Wherein, B is the strain matrix of the triangular mesh unit, and B=LN, L is the differential operator in the plane problem, and N is the 2×6 shape function matrix; K e is the stiffness matrix of the triangular mesh element; D is the elastic coefficient matrix of the triangular mesh unit, expressed as: Among them, G xy is the shear modulus between the two material directions x and y; E x 、E y Respectively refers to the elastic modulus of the material along the x-axis and y-axis; μ xy Refers to the Poisson's ratio in the y-axis direction when tension is applied in the x-axis direction, μ yx Similarly; and the elastic modulus and Poisson's ratio satisfy the following relationship: And x μ yx =E y μ xy (4) The stress conversion between the coordinate axis and the material axis is as follows: Among them, σ x , σ y are the normal stresses of the unit in the x and y directions, τ xy is the shear stress of the unit in the xy plane; The strain conversion relationship between the main direction of the material and the coordinate axis direction is: Among them, ε x , ε y , γ xy is the normal strain and shear strain of the unit in the x and y directions; Substituting equation (7) back into equation (5) yields: in, Let D0 = T1(θ)DA -1 Substituting T1(-θ)A into equation (2), the deformation energy of the triangular mesh unit is obtained as follows: The total deformation energy of the wound surface mesh model is the sum of the deformation energies of each triangular mesh unit, that is:

2. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 1 is characterized in that: The wound surface flattening shape optimization model is established according to the equivalent node load, with the total deformation potential energy of the unit material model as the objective function, and the optimized displacement of each unit node on the wound surface is solved, including: In order to describe the loading of the equivalent node load on the grid unit nodes before and after the expansion, the general form of the potential energy functional of the triangular grid unit is established as: Among them, P e is the unit equivalent node load matrix, and P e =K e δ e ; Introduce new unit node displacements into the wound surface mesh model: δ=[u1 v1 … u i v i … u n v n ] T (12) Then the overall deformation energy of the model after expansion is: Among them, G is the transformation matrix, which realizes the unit stiffness matrix K e The same order integration with the global stiffness matrix K is expressed as: Where E is the 2nd-order identity matrix; According to the variational principle, the condition for the functional to take an extreme value is that its first variation is 0, so we get: Kδ=λP (14) in, λ is the adjustment amplitude of displacement, which is used to adjust the number of model optimizations. Substituting an appropriate λ into equation (14) yields the optimal displacement δ of each grid unit node on the wound surface.

3. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 1 is characterized in that: The wound surface flattening boundary obtained through fitting and optimization includes: The contour points on the optimized wound surface unfolding shape are the type value points of the unfolded skin contour curve, which are back-calculated through the non-uniform rational B-spline control points; the first and last type value points are determined by the method of heavy nodes, and the remaining contour points are used as segmentation points in the fitting curve. The shape of the curve between adjacent type value points is described in the form of parametric equations; the skin contour curve formed by the optimized displacement of each grid unit node is defined as the original curve S m , where m is the scale and the node vector is denoted as U = [0,0,0,0,u4,…,u m+2 ,1,1,1,1], the adjacent type value points are used as the minimum section in the B-spline curve, and the distance between adjacent type value points determines the distribution of the node vector. The standard accumulated chord length method is used to determine the next node vector. The calculation method is: Among them, q i Represented as data in a node vector; According to the NURBS basis function of each segment of the curve, the relationship between the control points and the shape value points is: Among them, u i+3 is the starting node, i∈[1,m-1]; C i (i=0,1,2,…,m+2) is the corresponding i-th control point, V i is the data point, B represents the spline basis function; since the skin wound contour is a typical closed curve, the first and last points are the same, that is, V0 = V m , a total of m-1 equations; the first and last three control points overlap in sequence, and there is C m =C0,C m+1 =C1,C m+2 =C2, the number of unknowns also becomes m-1; the number of equations is equal to the number of unknowns, so all the control points of the skin contour curve are obtained by solving the following formula, and the specific shape fitting of the curve is realized by using the obtained control points to obtain the optimized wound flattening boundary; Where P (i＝0,...,m-1) It is the type value point.

4. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 1, characterized in that: The topological relationship between points, edges and surfaces of the wound surface mesh model is established based on the wound surface data, including: Define point table, patch vertex table, edge table, face table and topological hash function; The geometric information of the points contained in the wound surface data is stored in the point table; The face vertex table is required to not only retrieve the coordinate information of the corresponding point, but also record the topological relationship between the point and the surrounding point-edge-face; adjacent points, neighboring edges, adjacent edges and neighboring triangles are introduced into the face vertex table to enrich the topological association of the vertices; Adding two vertex indices, neighboring points and adjacent triangle indices constituting the edge to the edge table; Add the normal vector, three vertex indices, three edge indices and adjacent triangle indices of the face patch to the face table; The establishment of the topological hash function utilizes the hash function f v The key value K of the point is converted into the corresponding hash table address to access the point index in the face vertex table. The conversion formula is: f v =KmodM Among them, x, y, z are the coordinate components of the vertex in each direction, M is the length of the hash table, which is a prime number; similarly, using the hash function f e The key value and hash address of the edge are calculated and stored in the edge table. The information of the triangular facets in the wound surface data is read in sequence, and the coordinate information of the normal vector and the vertex is stored in the face table to achieve rapid topological reconstruction of the wound surface mesh model.

5. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 1, characterized in that: The step of dividing the wound surface mesh model into mesh hierarchies based on the hierarchical value and determining the initial expansion base point comprises: The concept of the level value is derived from the description of the progressive change of the level of the surface, starting from the outer boundary of the surface and gradually approaching inward; each point is assigned a level value vl, and points with the same level value are located at the same level. The corresponding level contour on the mesh surface is extracted according to the level value; wherein, the edges of two adjacent triangles are internal edges, and only one adjacent triangle edge is a boundary edge, the two adjacent points of the boundary edge are boundary points, and the triangle containing the boundary edge is a boundary triangle, and the level value is used as a basis to realize the assignment and division of the mesh points; After traversing all points on the wound mesh surface, the maximum level value is obtained as vl max , re-assign the level values ​​of all vertices, the assignment formula is: vl = vl max -vl+1(19) When mapping a spatial surface to a given plane, the initial unfolding base point is the first point mapped to the plane, and the point where vl=1 is taken as the initial unfolding base point of the wound mesh model; when this point is not unique, the center position of the entire surface mesh is calculated, and the point closest to the center position is selected as the initial unfolding base point.

6. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 1, characterized in that: Starting from the initial expansion base point, adjacent grid points are expanded layer by layer from the inside to the outside to a given plane based on the equal deformation method to obtain the flattened shape of the wound surface, including: The concept of equal deformation originates from the assumption of uniform continuity of material deformation. For any three-dimensional surface, a plane with topological equivalent can be found in two-dimensional space. The mesh on the wound surface mesh model is regarded as a triangular mesh unit, and the triangular mesh unit of the wound surface mesh model before expansion is recorded as ΔA0B0C0, and the corresponding triangular mesh unit after expansion by the equal deformation method is recorded as ΔABC. The unit area change rate before and after expansion is defined as τ, then: Among them, S i0 S is the area of ​​the triangle formed by the three vertices except i before unfolding. i is the area of ​​the unfolded triangle formed by the three vertices except i; After the expansion, the area of ​​the triangular mesh unit can also be expressed by the corresponding three vertex coordinates (x i ,y i )express: Among them, vertex D is the adjacent point to be expanded on the plane, and the coordinates (x D ,y D )unknown; From the obtained unit area, we can infer that: S i =T i i=A,B,C(22) In order to calculate the optimal approximate solution of formula (22), the distance metric M is introduced: M(S,T)=(S A -T A ) 2 +(S B -T B ) 2 +(S C -T C ) 2 (23) Through M, the unit area obtained by equation (22) and equation (23) is closest, that is, the minimum distance between S and T under the distance metric M is solved. Since M is non-negative, the extreme point of M is solved to determine the minimum value of M. Combined with equation (20), the unknown number x is obtained D ,y D The expression is: Among them, vertex D0 is the adjacent point to be developed on the surface; so far, the plane coordinates (x D ,y D ) to realize the planar unfolding of the mesh surface.

7. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 3 is characterized in that: The method further includes using a multi-resolution smoothing algorithm to perform smoothing on the flattened boundary of the wound surface: Construct a low-pass filter group to filter out the high-frequency part of the original curve to obtain a smooth low-frequency curve. The algorithm is: n =[<Φ n |Φ n >] -1 [<Φ n |Φ m >]C m ; Among them, [<Φ n |Φ n >] -1 [<Φ n |Φ m >] A low-pass filter group constructed to achieve curve smoothing, C m , C n are the control points of the target curve before and after smoothing respectively; the control points are substituted into equation (17) to complete the drawing of the smoothed curve, namely:

8. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 5, characterized in that: The method further comprises: The developability of the surface is described by the following formula: Where A is the total area of ​​the mesh surface, A sum (P i ) represents point P i The total area of ​​all neighboring triangles, θ sum (P i ) is point P i The sum of the interior angles formed by all its neighboring triangles; when point P i When it is a boundary point, take θ sum (P i )=2π; D(F) is between 0 and 1. When D(F) is closer to 0, it means that the surface is less developable and more likely to deform when unfolded. The clustering threshold ε is set. When D(F) < ε, the surface is considered to have poor developability and no clustering is required. When D(F) ≥ ε, the surface is considered to have good developability and further clustering is required. The clustering threshold is based on the vertex Gaussian curvature value K in the wound surface mesh model. p The value range of is given; The clustering method uses K-means clustering, randomly selecting K points in the grid as the initial center of each class, integrating the remaining vertices into the nearest class according to their similarity with each cluster center, and iteratively updating the value of each cluster center until the expected clustering result is met.

9. The three-dimensional wound surface finite element development method for mapping skin material properties according to claim 6, characterized in that: The method further comprises: For the treatment of morphological distortion points, the concept of developability needs to be introduced in combination with the uniform continuity assumption of the equal deformation method. Then the developability function of the vertices in the triangular mesh unit is: Among them, θ i It refers to the inner angle formed by the i-th neighborhood triangle of a vertex and the vertex. The threshold η is set. When D(v)≤η, it means that the point has not undergone a mutation and is a regular point. When D(v)>η, it means that the point has undergone a mutation and is defined as a morphological distortion point, but does not include boundary points. When processing the expansion of morphological deformity points on the wound surface, the expansion mutation of the morphological deformity points is evenly spread to its adjacent points in the form of neighborhood angles; on the basis of the equal deformation method, the adjacent edge invariance principle is added to ensure that the length of each neighborhood edge of the morphological deformity point remains unchanged during the expansion, and the angle ratio formed by each adjacent edge of the morphological deformity point remains unchanged, and the morphological deformity point is placed at the center of the neighborhood polygon formed by its neighboring points.

Images