Digitalized intelligent mse automatic design method based on biomechanical simulation

Abstract

This invention provides a digital intelligent MSE automatic design method based on biomechanical simulation, belonging to the field of biomechanical simulation technology. The method includes: acquiring three-dimensional morphological data and physiological load data of the target anatomical structure and mapping them to obtain a load distribution matrix; performing tensor operations on the load distribution matrix and tissue mechanical response characteristics to obtain a stress prediction field; identifying high-stress concentration areas and low-stiffness weak areas based on the stress prediction field; extracting contour curves and calculating the normal gradient to determine the expansion constraint direction; minimizing the stress distribution variance and the total tissue deformation energy and determining the spatial configuration parameters of the MSE device; using the spatial configuration parameters as boundary conditions to iteratively calculate the displacement field and strain field; correcting the spatial configuration parameters based on the difference between the cumulative deformation and a preset physiological threshold and the iterative deformation energy change rate; and outputting the corrected configuration parameters as design instructions when the convergence condition is met.

Description

Technical Field

[0001] This invention relates to the field of biomechanical simulation technology, and in particular to a digital intelligent MSE automatic design method based on biomechanical simulation. Background Technology

[0002] In the field of orthodontics and maxillofacial surgery, the design of MSE devices has long relied on the experience and manual measurement of clinicians. Conventional techniques involve selecting fixed-size pre-made or semi-finished expanders based on the patient's two-dimensional cephalometric data and standardized anatomical parameters. Through trial fitting and intraoperative adjustments to adapt to individual differences, doctors typically determine the implantation position and direction of the expander screws based on static anatomical landmarks and assess bone resistance strength through subjective palpation and imaging observation.

[0003] Current conventional techniques lack quantitative analysis of biomechanical response mechanisms, considering only geometric matching of anatomical morphology and failing to incorporate the stress-strain distribution of bone and soft tissues under physiological loads. This makes it difficult to accurately predict stress concentration areas when applying expansion forces, easily leading to complications such as local cortical bone cracking or root resorption. Secondly, adjusting design parameters relies on empirical trial and error. When abnormal expansion resistance or asymmetrical displacement occurs after the initial implantation, doctors need to gradually correct the screw angle or force frequency through multiple follow-up visits, extending the treatment course and making it impossible to guarantee that each correction will align with the optimal mechanical path. This can easily lead to cumulative and uncontrolled tissue deformation, affecting the stability and predictability of the orthodontic outcome. Summary of the Invention

[0004] This invention provides a digital intelligent MSE automatic design method based on biomechanical simulation, which can at least solve some of the problems existing in the prior art.

[0005] A first aspect of this invention provides a digital intelligent MSE automatic design method based on biomechanical simulation, comprising:

[0006] The three-dimensional morphological data and physiological load data of the target anatomical structure are acquired. The physiological load data is mapped to the mechanical bearing path corresponding to the three-dimensional morphological data to obtain the load distribution matrix. Tensor operations are performed on the load distribution matrix and the pre-stored tissue mechanical response characteristics to obtain the stress prediction field.

[0007] Based on the stress prediction field, high stress concentration areas and low stiffness weak areas are identified. Boundary tracking is performed on the high stress concentration areas and contour curves are extracted. The normal gradient is calculated along the contour curves and the expansion constraint direction is determined. Based on the expansion constraint direction and the low stiffness weak areas, the stress distribution variance and the total structural deformation energy are minimized and the spatial configuration parameters of the MSE device are determined.

[0008] The spatial configuration parameters are used as boundary conditions to iteratively calculate the displacement field and strain field. The displacement field is integrated over time to obtain the cumulative deformation, and the strain field is integrated spatially to obtain the iterative deformation energy. Based on the difference between the cumulative deformation and a preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined, and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is determined whether the iterative deformation energy meets the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions.

[0009] In one alternative implementation,

[0010] The process involves acquiring three-dimensional morphological data and physiological load data of the target anatomical structure, mapping the physiological load data onto the mechanical bearing path corresponding to the three-dimensional morphological data to obtain a load distribution matrix, and performing tensor operations on the load distribution matrix and pre-stored tissue mechanical response features to obtain a stress prediction field, including:

[0011] The tomographic image sequence of the target anatomical structure is acquired by medical imaging scanning equipment. The tomographic image sequence is reconstructed in three dimensions to obtain three-dimensional morphological data and the bone region and soft tissue region are segmented and extracted. The stress data of the target anatomical structure during physiological activities is acquired to obtain physiological load data and the direction is decomposed to obtain normal load component and tangential load component.

[0012] Feature points are extracted based on the boundary of the skeletal region and a skeletal topological connection relationship is constructed. The stiffness distribution field is calculated based on the thickness distribution of the soft tissue region. The mechanical bearing path is determined based on the skeletal topological connection relationship and the stiffness distribution field.

[0013] The normal load component is projected along the normal direction of the mechanical bearing path and accumulated to obtain the normal load distribution. The tangential load component is projected along the tangential direction of the mechanical bearing path and accumulated to obtain the tangential load distribution. A load distribution matrix is ​​constructed based on the normal load distribution and the tangential load distribution.

[0014] The elastic modulus and Poisson's ratio are extracted from the mechanical response characteristics of the tissue. A stiffness tensor is constructed based on the elastic modulus and Poisson's ratio. The load distribution matrix is ​​transposed and tensor contraction is performed with the stiffness tensor to obtain the initial stress tensor. The principal stress components and shear stress components of the initial stress tensor are calculated and a stress prediction field is constructed.

[0015] In one alternative implementation,

[0016] Based on the stress prediction field, high stress concentration regions and low stiffness weak regions are identified. Boundary tracing and contour curve extraction are performed on the high stress concentration regions. The normal gradient is calculated along the contour curve, and the expansion constraint direction is determined, including:

[0017] Morphological expansion and corrosion operations are performed on the stress prediction field to extract local extreme value regions. The stress mean of the local extreme value regions is calculated and compared with the global stress mean to determine high stress concentration regions. Stiffness feature values ​​of each voxel are extracted from the three-dimensional morphological data and a stiffness field is constructed. The regions in the stiffness field with stiffness feature values ​​lower than a preset stiffness threshold and the regions in the stress prediction field with stress values ​​lower than a preset stress threshold are spatially intersected to obtain low stiffness weak regions.

[0018] Voxelization boundary extraction is performed on the high stress concentration region to obtain a boundary voxel set. A three-dimensional point cloud is constructed based on the boundary voxel set and surface fitting is performed to obtain a boundary surface. Isoparametric sampling is performed on the boundary surface to obtain a sampling point sequence and a contour curve is constructed along the connection direction.

[0019] At each sampling point on the contour curve, the unit normal vector of the boundary surface is calculated. The directional derivative of the stress prediction field is calculated along the direction of the unit normal vector to obtain the normal stress gradient. The normal stress gradient is integrated along the contour curve, and the normal vector with the largest integral value is extracted as the principal expansion direction. The curvature tensor of the contour curve in three-dimensional space is calculated, and the principal direction corresponding to the minimum curvature is extracted as the secondary expansion direction. The expansion constraint direction is constructed based on the principal expansion direction and the secondary expansion direction.

[0020] In one alternative implementation,

[0021] Based on the expansion constraint direction and the low-stiffness weak region, minimizing the stress distribution variance and the total structural deformation energy, and determining the spatial configuration parameters of the MSE device, includes:

[0022] Extract the unit vector of the expansion constraint direction and construct the direction matrix. Determine the centroid coordinates of the low stiffness weak region and the high stress concentration region respectively and calculate the connection vector. Calculate the dot product of the connection vector and the direction matrix to obtain the direction correlation. Based on the direction correlation, determine the spatial position type of the low stiffness weak region relative to the expansion constraint direction.

[0023] Based on the spatial location type, candidate configuration regions are delineated in the three-dimensional morphological data and a set of sampling points is generated. The shortest distance between each sampling point in the set of sampling points and the boundary of the high stress concentration region is calculated to obtain the distance field. Combined with the expansion constraint direction, a feasible expansion path is determined. The path integral of the stress prediction field is performed along the feasible expansion path to obtain the cumulative path stress value. Combined with the length of the feasible expansion path, the average path stress is determined. Based on the average path stress, the change in stress distribution variance is determined.

[0024] The strain energy density is calculated by extracting the elastic modulus distribution of the three-dimensional morphological data. The strain energy density is integrated within the volume covered by the feasible expansion path to obtain the total tissue deformation energy. The comprehensive evaluation index is calculated by combining the change in the stress distribution variance and the directional correlation. The sampling point with the minimum comprehensive evaluation index is identified in the candidate configuration area, and its coordinates and direction vectors are extracted as spatial configuration parameters.

[0025] In one alternative implementation,

[0026] Extracting the elastic modulus distribution from the three-dimensional morphological data to calculate the strain energy density, and integrating the strain energy density over the volume covered by the feasible expansion path to obtain the total tissue deformation energy, including:

[0027] Multi-scale tissue layering structure is extracted from the three-dimensional morphological data and fiber orientation tensor reconstruction is performed on each layer. An anisotropic elastic tensor field is constructed based on the fiber orientation tensor and strain tensor components are calculated by combining the stress prediction field. Bilinear tensor contraction operation is performed on the strain tensor components and the anisotropic elastic tensor field to obtain the direction-dependent strain energy density field.

[0028] The geodesic distance field of the feasible expansion path is calculated and the path influence domain is constructed. Within the path influence domain, the direction-dependent strain energy density field is non-uniformly weighted and an adaptive integral grid is constructed. At each grid node, the local gradient of the direction-dependent strain energy density field is calculated and divergence is performed to obtain the energy flux density.

[0029] The boundary energy flux is obtained by performing an inner product operation between the energy flux density and the boundary normal vector of the path influence domain. The boundary energy flux is then subjected to surface integral along the boundary of the path influence domain and combined with the volume integral within the adaptive integral grid to obtain the path-related deformation energy. The path-related deformation energy corresponding to all feasible expansion paths is then nonlinearly aggregated to obtain the total tissue deformation energy.

[0030] In one alternative implementation,

[0031] The displacement field and strain field are obtained by iterative calculation using the spatial configuration parameters as boundary conditions. The cumulative deformation is obtained by time integration of the displacement field, and the iterative deformation energy is obtained by spatial integration of the strain field, including:

[0032] The position coordinates and expansion direction parameters of the MSE device are extracted from the spatial configuration parameters. Based on the position coordinates, the position of the point of action is calibrated in the three-dimensional morphological data and the force boundary conditions are constructed in combination with the expansion direction parameters. Fixed constraint regions are identified from the three-dimensional morphological data. Zero displacement constraints are applied to the fixed constraint regions to obtain displacement boundary conditions and combined with the force boundary conditions to obtain complete boundary conditions.

[0033] The three-dimensional morphological data is discretized using the finite element method to construct an overall stiffness matrix. The complete boundary conditions are applied to the overall stiffness matrix to obtain the displacement vectors of each node. Based on the displacement vectors, shape function interpolation is performed to obtain a continuous displacement field. The spatial partial derivatives of the displacement vectors of each voxel in the continuous displacement field are calculated to obtain the displacement gradient tensor. The displacement gradient tensor is symmetric to obtain the strain tensor and arranged to obtain the strain field. The displacement increment vectors of the continuous displacement field between adjacent time steps are calculated and accumulated along the time axis using the trapezoidal integral method to obtain the cumulative deformation.

[0034] Anisotropic elastic tensors are extracted from the mechanical response characteristics of the tissue and bilinearly contracted with the strain tensors of each voxel in the strain field to obtain the voxel strain energy density. The voxel strain energy density is then volume integraled in the spatial domain to obtain the iterative deformation energy.

[0035] In one alternative implementation,

[0036] Based on the difference between the cumulative deformation amount and the preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined, and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is then determined whether the iterative deformation energy satisfies the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions, including:

[0037] The deformation deviation is obtained by calculating the difference between the cumulative deformation and the preset physiological threshold. The deformation deviation is then judged by sign and marked as exceeding the limit. The difference between the iterative deformation energy and the iterative deformation energy corresponding to the previous iteration is calculated, and the rate of change of iterative deformation energy is calculated based on the difference.

[0038] Logarithmic transformation is performed on the absolute value of the deformation deviation to obtain the deviation scale, and reciprocal transformation is performed on the absolute value of the rate of change of the iterative deformation energy to obtain the convergence rate. Hyperbolic tangent mapping is performed on the deviation scale and the convergence rate to obtain the adaptive step size factor, and the adjustment step size is calculated by combining the preset benchmark step size.

[0039] The gradient vector is obtained by calculating the partial derivatives of the iterative deformation energy with respect to each parameter component in the spatial configuration parameters. The gradient vector is normalized by the second norm to determine the unit gradient direction. The correction direction coefficient is determined based on the over-limit state, and the correction increment vector is calculated by combining the unit gradient direction and the adjustment step size. The correction increment vector is superimposed on the spatial configuration parameters to obtain the correction configuration parameters. It is determined whether the absolute value of the rate of change of the iterative deformation energy is less than a preset convergence threshold. If it is less, the correction configuration parameters are output as design instructions.

[0040] A second aspect of this invention provides a digital intelligent MSE automatic design system based on biomechanical simulation, comprising:

[0041] The mechanical analysis unit is used to acquire three-dimensional morphological data and physiological load data of the target anatomical structure, map the physiological load data to the mechanical bearing path corresponding to the three-dimensional morphological data to obtain a load distribution matrix, and perform tensor operations on the load distribution matrix and pre-stored tissue mechanical response characteristics to obtain a stress prediction field.

[0042] The parameter optimization unit is used to identify high stress concentration areas and low stiffness weak areas based on the stress prediction field, perform boundary tracking and extract the contour curve of the high stress concentration area, calculate the normal gradient along the contour curve and determine the expansion constraint direction, and minimize the stress distribution variance and the total structural deformation energy according to the expansion constraint direction and the low stiffness weak area, and determine the spatial configuration parameters of the MSE device.

[0043] The iteration involves a unit used to iteratively calculate the displacement field and strain field by using the spatial configuration parameters as boundary conditions. The displacement field is integrated over time to obtain the cumulative deformation, and the strain field is integrated spatially to obtain the iterative deformation energy. Based on the difference between the cumulative deformation and a preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is then determined whether the iterative deformation energy meets the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions.

[0044] A third aspect of the present invention provides an electronic device, comprising:

[0045] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0046] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0047] In this invention, the fusion of tensor operations and load distribution matrix significantly improves the calculation accuracy of the stress prediction field. It can automatically identify high stress concentration areas and low stiffness weak areas, and determine the expansion constraint direction along the normal gradient of the contour curve. Thus, in the process of minimizing the stress distribution variance and the total tissue deformation energy, it achieves precise optimization of the spatial configuration parameters of the MSE device, effectively reducing the risk of postoperative stress concentration and the probability of tissue damage. Using the spatial configuration parameters as boundary conditions for iterative calculation, it can dynamically generate displacement and strain fields, and obtain the cumulative deformation and iterative deformation energy through time integration and spatial integration, respectively. Based on the difference between the cumulative deformation and the physiological threshold, combined with the rate of change of iterative deformation energy in continuous iterations, it adaptively adjusts the step size and corrects the configuration parameters along the gradient direction, which greatly improves the convergence speed and stability of the design process, avoids the repetitive calculations and resource waste caused by manual trial and error, and can automatically output corrected configuration parameters as design instructions by judging whether the iterative deformation energy meets the preset convergence conditions. It realizes full-process automation from data input to parameter output, significantly shortens the design cycle, reduces human error, and ensures that the final parameters achieve optimal balance in biomechanical performance, providing reliable technical support for personalized medicine. Attached Figure Description

[0048] Figure 1 This is a flowchart illustrating the digital intelligent MSE automatic design method based on biomechanical simulation, as described in an embodiment of the present invention.

[0049] Figure 2 This is a flowchart illustrating the deformation energy calculation process of an MSE device based on a digital intelligent MSE automatic design method using biomechanical simulation, as described in an embodiment of the present invention. Detailed Implementation

[0050] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0051] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0052] Figure 1 This is a flowchart illustrating the digital intelligent MSE automatic design method based on biomechanical simulation, as described in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0053] The three-dimensional morphological data and physiological load data of the target anatomical structure are acquired. The physiological load data is mapped to the mechanical bearing path corresponding to the three-dimensional morphological data to obtain the load distribution matrix. Tensor operations are performed on the load distribution matrix and the pre-stored tissue mechanical response characteristics to obtain the stress prediction field.

[0054] Based on the stress prediction field, high stress concentration areas and low stiffness weak areas are identified. Boundary tracking is performed on the high stress concentration areas and contour curves are extracted. The normal gradient is calculated along the contour curves and the expansion constraint direction is determined. Based on the expansion constraint direction and the low stiffness weak areas, the stress distribution variance and the total structural deformation energy are minimized and the spatial configuration parameters of the MSE device are determined.

[0055] The spatial configuration parameters are used as boundary conditions to iteratively calculate the displacement field and strain field. The displacement field is integrated over time to obtain the cumulative deformation, and the strain field is integrated spatially to obtain the iterative deformation energy. Based on the difference between the cumulative deformation and a preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined, and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is determined whether the iterative deformation energy meets the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions.

[0056] In one alternative implementation,

[0057] The process involves acquiring three-dimensional morphological data and physiological load data of the target anatomical structure, mapping the physiological load data onto the mechanical bearing path corresponding to the three-dimensional morphological data to obtain a load distribution matrix, and performing tensor operations on the load distribution matrix and pre-stored tissue mechanical response features to obtain a stress prediction field, including:

[0058] The tomographic image sequence of the target anatomical structure is acquired by medical imaging scanning equipment. The tomographic image sequence is reconstructed in three dimensions to obtain three-dimensional morphological data and the bone region and soft tissue region are segmented and extracted. The stress data of the target anatomical structure during physiological activities is acquired to obtain physiological load data and the direction is decomposed to obtain normal load component and tangential load component.

[0059] Feature points are extracted based on the boundary of the skeletal region and a skeletal topological connection relationship is constructed. The stiffness distribution field is calculated based on the thickness distribution of the soft tissue region. The mechanical bearing path is determined based on the skeletal topological connection relationship and the stiffness distribution field.

[0060] The normal load component is projected along the normal direction of the mechanical bearing path and accumulated to obtain the normal load distribution. The tangential load component is projected along the tangential direction of the mechanical bearing path and accumulated to obtain the tangential load distribution. A load distribution matrix is ​​constructed based on the normal load distribution and the tangential load distribution.

[0061] The elastic modulus and Poisson's ratio are extracted from the mechanical response characteristics of the tissue. A stiffness tensor is constructed based on the elastic modulus and Poisson's ratio. The load distribution matrix is ​​transposed and tensor contraction is performed with the stiffness tensor to obtain the initial stress tensor. The principal stress components and shear stress components of the initial stress tensor are calculated and a stress prediction field is constructed.

[0062] Medical imaging scanning equipment can be cone-beam computed tomography (CBCT), multi-slice spiral CT, or magnetic resonance imaging (MRI). Taking CBCT as an example, the equipment rotates 360 degrees around the patient's head, acquiring one frame of two-dimensional projection image every 0.3 to 0.5 degrees, obtaining 720 to 1200 tomographic images in a complete scan cycle. The slice thickness of these tomographic images is typically 0.2 to 0.5 millimeters, with a pixel resolution of 0.1 to 0.3 millimeters, which can clearly display the fine structure of the craniofacial bones and the morphological features of the surrounding soft tissues.

[0063] After obtaining the tomographic image sequence, a method combining grayscale thresholding segmentation and region growing was used for 3D reconstruction. First, grayscale normalization was performed on each frame of the tomographic image, mapping pixel grayscale values ​​to a standard range of 0 to 255. Then, the grayscale threshold range for skeletal tissue was set to 150 to 255, and for soft tissue, it was set to 30 to 150, thus initially distinguishing between skeletal and soft tissue regions through thresholding segmentation. To eliminate noise interference and isolated pixels, morphological opening and closing operations were performed on the segmentation results. The structuring element size for opening was 3×3 pixels, and for closing it was 5×5 pixels. Next, a region growing algorithm was used for fine extraction of the skeletal region. The center of the bone was selected as the seed point, and the growth criterion was set as the grayscale difference between adjacent pixels being less than 20 and within the bone threshold range. Iterative growth was performed until the region boundary stabilized. For the soft tissue region, the same region growing method was used, but the growth criterion was adjusted to the grayscale difference between adjacent pixels being less than 15 and within the soft tissue threshold range. After segmenting all tomographic images, three-dimensional voxels are stacked according to the spatial relationship of the images. The moving cube algorithm is used to extract the bone surface and soft tissue surface, and a triangular mesh model is generated as three-dimensional morphological data.

[0064] Acquiring physiological load data requires simulating the stress state of the target anatomical structure during actual physiological activities. For the craniofacial skeletal system, the main physiological loads originate from the contractile force of the masticatory muscles, occlusal contact force, and the constraint reaction force of the temporomandibular joint. Electromyography (EMG) is used to record the EMG signals of the masticatory muscles under different occlusal states. Based on the linear relationship between the EMG signal amplitude and muscle contractile force, the contractile forces of the masseter, temporalis, medial pterygoid, and lateral pterygoid muscles are calculated. Typical value ranges are: masseter 200-400 N, temporalis 150-300 N, medial pterygoid 100-200 N, and lateral pterygoid 50-100 N. Occlusal contact force is measured using a pressure-sensitive membrane. A 0.1 mm thick pressure membrane is placed between the upper and lower teeth, and the pressure distribution at each contact point during occlusion is recorded. The pressure range at a single contact point is 10-50 N. The constraint reaction force of the temporomandibular joint is obtained through inverse dynamics calculation. Based on the movement trajectory of the mandible and the muscle contraction force, the support reaction force at the joint is solved, with typical values ​​ranging from 50 Newtons to 150 Newtons.

[0065] After obtaining the force data, directional decomposition is required for subsequent load mapping. A three-dimensional Cartesian coordinate system is established with the centroid of the target anatomical structure as the origin, the X-axis pointing to the left, the Y-axis pointing forward, and the Z-axis pointing upward. For each load vector, its components in the three coordinate axes are calculated. The normal load component is defined as the projection of the load vector onto the local surface normal direction, calculated by taking the dot product of the load vector and the surface normal vector. The tangential load component is defined as the projection of the load vector onto the local surface tangent plane, calculated by subtracting the normal load component from the load vector; the remaining part is the tangential load component. For the contractile force of the masticatory muscles, the direction of action is determined by the spatial position of the muscle's origin and insertion points, represented by a vector connecting the origin and insertion points. For the occlusal contact force, its direction of action is perpendicular to the occlusal surface; therefore, the normal load component is dominant, while the tangential load component is smaller. For the constraint reaction force of the temporomandibular joint, its direction of action is determined by the morphological characteristics of the glenoid fossa, typically forming an angle of 15 to 30 degrees with the normal of the articular surface.

[0066] The mechanical load-bearing path is determined by comprehensively considering the topological structure of the skeleton and the stiffness distribution of soft tissue. For the skeletal region, feature points are extracted from the bone surface. The selection principle for feature points is locations with significant curvature changes, including anatomical landmarks such as bony prominences, ridges, and edges of bone foramina. A principal curvature calculation method is used to calculate the principal curvature value for each vertex on the bone surface. Vertices with an absolute principal curvature value greater than 0.05 are selected as candidate feature points. Then, a non-maximum suppression algorithm is used to filter out local curvature extrema as the final feature points. The number of feature points is typically 200 to 500. Based on the aforementioned feature points, a skeletal topological connection relationship is constructed. The Delaunay triangulation algorithm is used to establish connecting edges between feature points, forming a skeletal topological network. The skeletal topological network reflects the main force transmission paths within the bone, and the weight of the connecting edge is determined by the geodesic distance between two feature points and the average bone density.

[0067] For soft tissue regions, the stiffness distribution field is calculated based on the thickness distribution of the soft tissue. The soft tissue thickness is obtained by calculating the shortest distance from the soft tissue surface to the bone surface; the thickness typically ranges from 1 mm to 10 mm. The local stiffness of soft tissue is inversely proportional to its thickness; regions with greater thickness have lower stiffness, and regions with less thickness have higher stiffness. In specific calculations, a reference stiffness value of 1 MPa is set, corresponding to a reference thickness of 5 mm. When the actual thickness deviates from the reference thickness, the stiffness value is adjusted inversely proportionally. For example, the stiffness of a soft tissue region with a thickness of 2 mm is approximately 2.5 MPa, and the stiffness of a soft tissue region with a thickness of 8 mm is approximately 0.625 MPa. The entire soft tissue region is divided into several voxel units, each with dimensions of 0.5 mm × 0.5 mm × 0.5 mm. A corresponding stiffness value is assigned to each voxel unit, forming the stiffness distribution field.

[0068] Starting from the load application point, a path search is performed along the edge with the minimum weight in the skeletal topology network, while considering the influence of soft tissue stiffness distribution. When the path passes through soft tissue regions, paths with higher stiffness are prioritized to avoid weak areas with excessively low stiffness. Dijkstra's shortest path algorithm is used for path planning, combining the edge weights of the skeletal topology network with the reciprocal of the soft tissue stiffness to construct a comprehensive weight function. The final mechanical bearing path is the optimal force transmission path from the load application point to the constraint support point, where each node has explicit spatial coordinates and local stiffness information.

[0069] The normal and tangential load components are mapped onto the mechanical bearing path to construct a load distribution matrix. For the normal load component, the normal projection distance from the load application point to each node on the mechanical bearing path is calculated. A Gaussian attenuation function is used to describe the spatial distribution characteristics of the load, with an attenuation coefficient set to 0.2, meaning that for every 1 mm increase in distance from the load application point, the load amplitude attenuates to 81.87% of its original value. The normal load component is distributed to each node on the path according to the Gaussian attenuation law, and the contribution of all load application points to the same node is accumulated to obtain the normal load distribution. For the tangential load component, a similar method is used for tangential projection and accumulation, but the attenuation coefficient is adjusted to 0.15 to reflect the relatively wide transmission range of the tangential load. The normal and tangential load distributions are arranged in node order to form a load distribution matrix. The number of rows in the matrix is ​​equal to the number of nodes on the mechanical bearing path, and the number of columns is 2. The first column stores the normal load distribution, and the second column stores the tangential load distribution.

[0070] The tissue mechanical response characteristics encompass the material properties of each component of the target anatomical structure. For skeletal tissue, the elastic modulus is determined based on bone density; the elastic modulus of cortical bone is 15-20 GPa, and that of cancellous bone is 0.5-2 GPa, with a Poisson's ratio consistently set at 0.3. For soft tissue, the elastic modulus is determined based on tissue type; the elastic modulus of gingival mucosa is 0.5-1 MPa, and that of periodontal ligament is 0.05-0.2 MPa, with a Poisson's ratio ranging from 0.45 to 0.49, approaching that of incompressible materials. A stiffness tensor is constructed based on these material properties. The stiffness tensor is a fourth-order tensor, which, under the isotropic assumption, can be simplified to a function of two independent parameters. The specific form of the stiffness tensor is determined by the elastic modulus and Poisson's ratio through constitutive relations, reflecting the linear relationship between stress and strain when the material is subjected to force.

[0071] The stress prediction field is calculated using tensor contraction. The load distribution matrix is ​​transposed to obtain a transpose matrix with two rows and columns equal to the number of nodes. Then, a tensor contraction operation is performed between the transpose matrix and the stiffness tensor, shrinking each column of the transpose matrix to obtain the initial stress tensor at the corresponding node. The initial stress tensor is a second-order symmetric tensor containing six independent components: three normal stress components and three shear stress components. Eigenvalue decomposition of the initial stress tensor yields three principal stress components, representing the maximum, intermediate, and minimum principal stresses, with the directions of the principal stresses corresponding to the directions of the eigenvectors. The shear stress components are directly obtained from the off-diagonal elements of the stress tensor. The principal stress components and shear stress components of all nodes are interpolated according to their spatial locations to generate a continuous stress prediction field. This field describes the stress distribution of the target anatomical structure under physiological loads, providing a mechanical basis for subsequent identification of high-stress concentration areas and the design of the MSE device.

[0072] In one alternative implementation,

[0073] Based on the stress prediction field, high stress concentration regions and low stiffness weak regions are identified. Boundary tracing and contour curve extraction are performed on the high stress concentration regions. The normal gradient is calculated along the contour curve, and the expansion constraint direction is determined, including:

[0074] Morphological expansion and corrosion operations are performed on the stress prediction field to extract local extreme value regions. The stress mean of the local extreme value regions is calculated and compared with the global stress mean to determine high stress concentration regions. Stiffness feature values ​​of each voxel are extracted from the three-dimensional morphological data and a stiffness field is constructed. The regions in the stiffness field with stiffness feature values ​​lower than a preset stiffness threshold and the regions in the stress prediction field with stress values ​​lower than a preset stress threshold are spatially intersected to obtain low stiffness weak regions.

[0075] Voxelization boundary extraction is performed on the high stress concentration region to obtain a boundary voxel set. A three-dimensional point cloud is constructed based on the boundary voxel set and surface fitting is performed to obtain a boundary surface. Isoparametric sampling is performed on the boundary surface to obtain a sampling point sequence and a contour curve is constructed along the connection direction.

[0076] At each sampling point on the contour curve, the unit normal vector of the boundary surface is calculated. The directional derivative of the stress prediction field is calculated along the direction of the unit normal vector to obtain the normal stress gradient. The normal stress gradient is integrated along the contour curve, and the normal vector with the largest integral value is extracted as the principal expansion direction. The curvature tensor of the contour curve in three-dimensional space is calculated, and the principal direction corresponding to the minimum curvature is extracted as the secondary expansion direction. The expansion constraint direction is constructed based on the principal expansion direction and the secondary expansion direction.

[0077] After obtaining the stress prediction field, it is necessary to accurately identify the regions that have a critical impact on the MSE device design. Morphological expansion and erosion operations are performed on the stress prediction field, controlling the spatial range of the operation by setting the size and shape of the structural elements. Expansion expands high-stress regions outward, while erosion shrinks the boundaries; the combination of these two operations effectively filters out noise and preserves the true stress concentration characteristics. After completing the morphological operations, the processed stress field is traversed, extracting all sets of voxels that satisfy the local extremum condition; the stress values ​​of these voxels reach their maximum within their neighborhood. For each extracted local extremum region, the stress values ​​of all voxels within that region are statistically analyzed, and the arithmetic mean is calculated as the representative stress of that region. Simultaneously, the global stress mean of the entire stress prediction field is calculated. The stress mean of each local extremum region is numerically compared with the global stress mean. When the local stress mean exceeds 1.5 to 2 times the global stress mean, the region is marked as a high-stress concentration region.

[0078] To identify low-stiffness weak areas, tissue density, elastic modulus, or other mechanical property parameters corresponding to each voxel location are extracted from the 3D morphological data. Stiffness eigenvalues ​​of each voxel are calculated based on these parameters; these eigenvalues ​​reflect the tissue's resistance to deformation at that location. All voxel stiffness eigenvalues ​​are organized into a 3D matrix according to their spatial location to construct a complete stiffness field. A stiffness threshold is set, which can be determined based on the normal tissue stiffness range of the target anatomical structure, typically 60% to 70% of the average normal stiffness. The stiffness field is iterated, and all voxels with stiffness eigenvalues ​​below the threshold are marked as insufficient stiffness areas. Simultaneously, a stress threshold is set in the stress prediction field, typically 0.3 to 0.5 times the global stress average; voxels with stress values ​​below this threshold are marked as low-stress areas. A 3D Boolean intersection operation is performed on the insufficient stiffness areas and the low-stress areas to extract the set of voxels that simultaneously satisfy both conditions. These voxels, which experience both low stress and low stiffness, constitute the low-stiffness weak areas requiring focused attention.

[0079] After identifying the high-stress concentration region, precise extraction of the spatial boundary is required to guide the positioning of the MSE device. Voxelized boundary extraction is performed on the high-stress concentration region. All voxels within the region are traversed, and the 26-neighborhood or 6-neighborhood of each voxel is examined. If at least one neighboring voxel does not belong to the high-stress concentration region, the voxel is marked as a boundary voxel. The 3D coordinates of all boundary voxels are collected to form a boundary voxel set. The center coordinates of each voxel in the boundary voxel set are used as spatial points to construct a 3D point cloud data structure. Surface reconstruction is performed on this point cloud using methods such as moving least squares, Poisson surface reconstruction, or radial basis function interpolation to fit a smooth boundary surface. This boundary surface is represented by parametric equations, has two parameter dimensions, and can fully describe the 3D geometry of the high-stress concentration region.

[0080] Isoparametric sampling is performed on the boundary surface, with sampling intervals set along the first and second parameter directions. The selection of the sampling interval needs to balance computational accuracy and efficiency, typically between 0.5 mm and 2 mm. A regular mesh is generated in the parameter space based on the sampling interval, and the mesh nodes are mapped to the three-dimensional physical space, resulting in a sequence of sampling points uniformly distributed on the boundary surface. Adjacent sampling points are connected sequentially according to the connectivity relationships in the parameter space, forming a closed or open spatial curve, which is the contour curve of the high-stress concentration region. The contour curve may contain multiple independent curve segments, each corresponding to a characteristic contour on the boundary surface.

[0081] The normal gradient along the contour curve is calculated to determine the preferred direction of tissue expansion. At each sampling point on the contour curve, the tangent plane at that point is calculated based on the parametric equations of the boundary surface. Two tangent vectors are obtained by taking the partial derivatives of the parametric equations, and the cross product of the two tangent vectors gives the normal vector at that point. The normal vector is normalized to obtain the unit normal vector, which is perpendicular to the boundary surface and points outward from the high stress concentration region. Starting from this sampling point, a series of neighboring points are selected in the stress prediction field along the direction of the unit normal vector. The stress values ​​at these points are calculated, and the rate of change of stress along the normal is calculated using the finite difference method to obtain the normal stress gradient at that sampling point. The magnitude of the normal stress gradient reflects the degree of stress change perpendicular to the boundary direction; a larger gradient value indicates a faster stress decay in that direction, which is a favorable direction for tissue expansion.

[0082] The normal stress gradient at all sampling points on the profile curve is integrated, with the integration path extending along the curve. The integration result reflects the cumulative effect of the normal stress gradient across the entire profile curve. The integral values ​​corresponding to the normal vectors at different sampling points are compared, and the normal vector with the largest integral value is extracted. This normal vector represents the direction with the most significant stress relief effect in an overall sense, and is determined as the main expansion direction. The main expansion direction is the primary direction in which the MSE device applies the expansion force; expansion along this direction can most effectively reduce the stress level in high stress concentration areas.

[0083] To comprehensively describe the spatial constraints of tissue expansion and determine the secondary expansion direction, the curvature characteristics of the contour curve in three-dimensional space are calculated. For each sampling point on the curve, the first and second derivatives of the curve are calculated based on the position of that point and its neighboring points. The curvature value at that point is obtained by calculating the cross product of the derivatives and the modulus. In three-dimensional space, the bending characteristics of the curve can be described by a curvature tensor, which contains information about the degree of curvature of the curve in different directions. Eigenvalue decomposition is performed on the curvature tensor to extract the eigenvector corresponding to the smallest eigenvalue. This eigenvector represents the direction of minimum curvature of the curve, i.e., the direction of the curve's flattest extension. This eigenvector is determined as the secondary expansion direction, which is usually orthogonal or approximately orthogonal to the primary expansion direction, together forming the spatial coordinate system of tissue expansion.

[0084] An expansion constraint direction is constructed based on the primary and secondary expansion directions. This constraint direction uses the primary expansion direction as the main axis and the secondary expansion direction as the auxiliary axis, forming a spatial reference frame for the force application of the MSE device. In practical applications, the weighting coefficient of the primary expansion direction is typically set to 0.7 to 0.85, and the weighting coefficient of the secondary expansion direction is set to 0.15 to 0.3. A comprehensive expansion constraint direction is obtained through weighted combination. This expansion constraint direction not only considers the stress release requirements but also takes into account the geometric characteristics of the anatomical structure, ensuring that the expansion action of the MSE device can effectively reduce stress concentration while adapting to the natural morphology of the tissue, avoiding the generation of new stress concentrations or tissue damage during the expansion process.

[0085] Based on the identified low-stiffness weak areas, when determining the specific installation location of the MSE device, it is necessary to avoid these areas or provide additional support in these areas. Due to insufficient tissue load-bearing capacity, directly applying expansion forces to these low-stiffness weak areas may lead to excessive local tissue deformation or even damage. By combining the spatial distribution information of low-stiffness weak areas with the direction of expansion constraints, the support point location of the MSE device can be optimized, ensuring that the support points fall in areas with higher stiffness. Simultaneously, the direction of the expansion force can effectively improve the stress state of the low-stiffness weak areas, promoting tissue reconstruction and mechanical property recovery in these areas.

[0086] In one alternative implementation,

[0087] Based on the expansion constraint direction and the low-stiffness weak region, minimizing the stress distribution variance and the total structural deformation energy, and determining the spatial configuration parameters of the MSE device, includes:

[0088] Extract the unit vector of the expansion constraint direction and construct the direction matrix. Determine the centroid coordinates of the low stiffness weak region and the high stress concentration region respectively and calculate the connection vector. Calculate the dot product of the connection vector and the direction matrix to obtain the direction correlation. Based on the direction correlation, determine the spatial position type of the low stiffness weak region relative to the expansion constraint direction.

[0089] Based on the spatial location type, candidate configuration regions are delineated in the three-dimensional morphological data and a set of sampling points is generated. The shortest distance between each sampling point in the set of sampling points and the boundary of the high stress concentration region is calculated to obtain the distance field. Combined with the expansion constraint direction, a feasible expansion path is determined. The path integral of the stress prediction field is performed along the feasible expansion path to obtain the cumulative path stress value. Combined with the length of the feasible expansion path, the average path stress is determined. Based on the average path stress, the change in stress distribution variance is determined.

[0090] The strain energy density is calculated by extracting the elastic modulus distribution of the three-dimensional morphological data. The strain energy density is integrated within the volume covered by the feasible expansion path to obtain the total tissue deformation energy. The comprehensive evaluation index is calculated by combining the change in the stress distribution variance and the directional correlation. The sampling point with the minimum comprehensive evaluation index is identified in the candidate configuration area, and its coordinates and direction vectors are extracted as spatial configuration parameters.

[0091] Determining the spatial configuration parameters of the MSE device requires comprehensive consideration of the spatial relationship between the expansion constraint direction and the low-stiffness weak region. The expansion constraint direction is vectorized, and the normal gradient calculated along the profile curve is normalized to a unit vector. Assuming that a high-stress concentration region is identified on the profile curve... At each of the key points, the normal gradient is normalized to obtain a unit vector. (in Arrange these unit vectors in spatial order to form a direction matrix. The dimension of this matrix is Each row represents a three-dimensional component of an expansion constraint direction.

[0092] For low-stiffness weak regions and high-stress concentration regions, the centroid coordinates are calculated separately. Low-stiffness weak regions are typically characterized by low bone density or insufficient cortical bone thickness. By thresholding the density distribution in the three-dimensional morphological data, the voxel set of this region can be identified. The centroid coordinates are obtained by arithmetically averaging the coordinates of all voxels within this set. Similarly, the centroid coordinates of high stress concentration regions The connection vector is obtained by averaging the coordinates of high-stress voxels exceeding a threshold in the stress prediction field. Defined as a vector pointing from the centroid of the high-stress-concentration region to the centroid of the low-stiffness-weak region, the calculation formula is: The vector is then normalized to obtain the unit connectivity vector. .

[0093] Directional dependence is calculated using the dot product operation. This involves connecting unit vectors. With direction matrix Each unit vector in Performing the dot product operation, we get scalar value The range of these scalar values ​​is within arrive Between these values, positive values ​​indicate that the connection vector and the expansion constraint are roughly in the same direction, while negative values ​​indicate the opposite direction, and the absolute values ​​are close to... This indicates a high degree of relevance. (Retrieve all) The maximum value is used as the overall directional correlation index. .according to The numerical range determines the spatial location type: when When it is determined to be "forward configuration type", it means that the low stiffness weak region is located on the extension path of the expansion constraint direction; when When it is determined to be "lateral configuration type", it means that diagonal support is required; when When the condition is determined to be "reverse configuration type", it means that the low stiffness weak area is located on the opposite side of the expansion constraint direction, requiring special mechanical transmission design.

[0094] Candidate configuration regions are delineated in 3D morphological data based on spatial location type. For forward configuration, candidate configuration regions are defined by connection vectors. With the central axis and radius as The candidate configuration region is a cylindrical region, the length of which is equal to the magnitude of the connection vector. For lateral configurations, the candidate configuration region is expanded into a cuboid region with the connection vector as its diagonal, the side length of which is adaptively determined according to the scale of the anatomical structure. For reverse configurations, the candidate configuration region is defined as an irregular region containing all possible force transmission paths between high-stress concentration regions and low-stiffness weak regions, obtained through morphological dilation. Within the candidate configuration region, a sampling point set is generated using uniform mesh sampling or a density-adaptive sampling strategy. The number of sampling points is determined based on the volume of the candidate region and the required accuracy, typically in... arrive Between one.

[0095] For each sampling point in the sampling point set The shortest distance between the sampling point and the boundary of the high stress concentration region is calculated. The boundary of the high stress concentration region is obtained by extracting isosurfaces from the stress prediction field, and the shortest distance from the sampling point to the boundary is calculated using either the fast traversal method or the signed distance field algorithm. Distance field It consists of the distance values ​​of all sampling points. When determining a feasible expansion path in conjunction with the expansion constraint direction, it starts from each sampling point. Start along the direction of the nearest expansion constraint. Perform ray tracing until the rays reach the boundary of the high-stress concentration area or extend beyond the anatomical structure. Feasible expansion path. Defined as a ray segment from the sampling point to the boundary, the feasibility of the path is determined by checking whether the path crosses prohibited areas (such as neurovascular bundles, tooth roots, etc.).

[0096] The stress prediction field is integraled along a feasible expansion path. For the path... Discretize it into For each line segment, the stress value of the stress prediction field is sampled at the midpoint of the stress prediction field. (in The cumulative path stress is calculated as follows: ,in For the first The length of each line segment. Path length. It is the sum of the lengths of all line segments, that is... The mean path stress is defined as... This indicator reflects the average stress level along the expansion path.

[0097] The change in stress distribution variance was calculated by comparing the stress distribution before and after the MSE device configuration. Assuming at the sampling point... After the MSE device is installed, the stress distribution along the feasible expansion path changes. The new stress distribution can be approximated by superimposing the original stress prediction field with the additional stress field generated by the MSE device. The magnitude of the additional stress field is related to the expansion force and direction of the MSE device, as well as the mechanical response characteristics of the tissue. The change in stress distribution variance... Calculate the difference between the variance of the configured stress field and the variance of the original predicted stress field. To simplify the calculation, an approximate formula based on average path stress can be used: ,in This is an empirical coefficient. This represents the maximum stress value in the stress prediction field. This formula shows that deploying an MSE device in high-stress regions can more effectively reduce the stress distribution variance.

[0098] The calculation of the total deformation energy of tissue is based on the theory of elasticity. The elastic modulus distribution is extracted from three-dimensional morphological data. This distribution reflects the stiffness differences among different tissue types; for example, the elastic modulus of cortical bone is approximately [missing value]. Cancellous bone is approximately Soft tissue is approximately Strain energy density Calculated under the assumption of online elasticity ,in The strain field can be approximated by dividing the predicted stress field by the elastic modulus. For sampling points... The corresponding feasible expansion path is defined, and the path coverage volume is defined. With the path as the central axis and a radius of The volume of the cylinder, where the radius represents the mechanical influence range of the MSE device, is typically taken as... arrive Millimeters. Total tissue deformation energy Through volume Integral calculation of internal strain energy density: The numerical implementation uses a voxel summation approximation: ,in The volume of a single voxel.

[0099] Comprehensive evaluation indicators The calculation requires balancing two objectives: optimizing stress distribution and minimizing tissue damage. A comprehensive evaluation index is defined as follows: ,in , , These are weighting coefficients, which respectively control the influence of stress distribution variance variation, tissue deformation energy, and orientation correlation. The selection of these weighting coefficients is based on clinical experience and optimization objectives; a typical value is [value to be filled in]. , , In this formula, A negative value indicates an improved stress distribution; therefore, multiplying by a positive weight makes... Decrease; A higher value indicates a higher tissue deformation energy, which should be avoided. The larger the value, the more consistent the configuration direction is with the direction of expansion constraints, which should be encouraged; hence the negative sign.

[0100] Traverse all sampling points within the candidate configuration area and calculate the comprehensive evaluation index for each sampling point. , recognition Sampling point location that reaches the minimum value This location is the optimal configuration position for the MSE device. Extract the three-dimensional coordinates of this location. The position component is used as a spatial configuration parameter. Simultaneously, the expansion constraint direction corresponding to this sampling point is extracted. The direction component serves as a spatial configuration parameter. The direction vector is typically represented as three components of a unit vector. Alternatively, Euler angles can be used. Furthermore, the length of the feasible extension path... Based on the mechanical response characteristics of the tissue, a preliminary estimate of the expansion required for the MSE device is made, serving as additional information for the spatial configuration parameters. These parameters together constitute a complete set of spatial configuration parameters, providing initial boundary conditions for subsequent iterative optimization.

[0101] In one alternative implementation,

[0102] Extracting the elastic modulus distribution from the three-dimensional morphological data to calculate the strain energy density, and integrating the strain energy density over the volume covered by the feasible expansion path to obtain the total tissue deformation energy, including:

[0103] Multi-scale tissue layering structure is extracted from the three-dimensional morphological data and fiber orientation tensor reconstruction is performed on each layer. An anisotropic elastic tensor field is constructed based on the fiber orientation tensor and strain tensor components are calculated by combining the stress prediction field. Bilinear tensor contraction operation is performed on the strain tensor components and the anisotropic elastic tensor field to obtain the direction-dependent strain energy density field.

[0104] The geodesic distance field of the feasible expansion path is calculated and the path influence domain is constructed. Within the path influence domain, the direction-dependent strain energy density field is non-uniformly weighted and an adaptive integral grid is constructed. At each grid node, the local gradient of the direction-dependent strain energy density field is calculated and divergence is performed to obtain the energy flux density.

[0105] The boundary energy flux is obtained by performing an inner product operation between the energy flux density and the boundary normal vector of the path influence domain. The boundary energy flux is then subjected to surface integral along the boundary of the path influence domain and combined with the volume integral within the adaptive integral grid to obtain the path-related deformation energy. The path-related deformation energy corresponding to all feasible expansion paths is then nonlinearly aggregated to obtain the total tissue deformation energy.

[0106] Tissue tomographic scan information was retrieved from stored 3D morphological data, and the spatial distribution of bone tissue, periodontal ligament, alveolar bone, and soft tissue was identified using a grayscale segmentation algorithm. For the bone tissue region, morphological dilation and erosion operations were used to extract the interface between the outer cortical bone and the inner cancellous bone. Gaussian filtering was further applied to the cortical bone region for noise reduction, and Hessian matrix eigenvalue analysis was used to identify the principal direction of trabeculae. Within the cancellous bone region, rays were projected along three orthogonal directions with a fixed step size. The number of trabeculae and the inter-ray length encountered on each ray were counted, and these statistical results were converted into a probability distribution of local fiber orientation. For the periodontal ligament tissue, by analyzing its attachment relationship with the root surface and the alveolar bone wall, the fiber orientation was set as a radial distribution pattern from the root surface to the alveolar bone wall.

[0107] Each identified tissue layer is discretized into three-dimensional voxels, and a local coordinate system is established at the center of each voxel. For voxels with fiber orientation characteristics, the fiber orientation vectors within the voxel and its neighborhood are collected to construct a second-order tensor. This tensor is expressed as a vector of fiber directions. The fiber orientation tensor is obtained by performing an outer product operation with itself, followed by a volume-weighted average of all fiber tensors in the neighborhood, resulting in a reconstructed tensor characterizing the statistical properties of fiber orientation at that location. For isotropic tissue regions, the fiber orientation tensor is set to a unit sphere tensor to ensure the continuity of subsequent mechanical calculations.

[0108] Based on the reconstructed fiber orientation tensor field, a fourth-order anisotropic elastic tensor is constructed at each voxel location. For transversely isotropic materials, the elastic tensor can be derived from the fiber-oriented elastic modulus. Vertical elastic modulus Poisson's ratio in the plane Poisson's ratio outside the surface and shear modulus This is defined by setting the fiber direction as the principal direction of the tensor and rotating the elastic tensor from the material's principal axes to the global coordinate system through coordinate transformation. For bone tissue, the elastic modulus parameter is determined based on the local bone density value using a power-law relationship, ensuring that the elastic modulus parameter in the cortical bone region is... Significantly higher than in the cancellous bone region. For the periodontal ligament, based on clinical trial data, its stiffness along the fiber direction is set to be approximately 5 to 8 times that in the vertical direction.

[0109] The stress tensor at each voxel location is extracted from the previously calculated stress prediction field. According to the tensor form of the generalized Hooke's law, this is expressed through the inverse tensor (compliance tensor) of the anisotropic elasticity tensor. The strain tensor is obtained by performing a contraction operation with the stress tensor. This operation, under the specified notation, is represented as the sum of the products of the compliance tensor components and the stress components. In the numerical implementation, the fourth-order compliance tensor is converted into a sixth-order square matrix using Voigt notation, and the stress tensor is converted into a six-dimensional column vector. The calculation is then performed using matrix multiplication. For each voxel, six independent components of the strain tensor are extracted, including three normal strain components. , , and three shear strain components , , .

[0110] When calculating the direction-dependent strain energy density, a bilinear tensor contraction operation is employed. This operation involves a combination of the strain tensor, elastic tensor, and stress tensor, specifically, the strain energy density is equal to half the dot product of the strain and stress tensors. Considering the anisotropic nature of the elastic tensor, the strain energy density exhibits significant differences in different spatial directions. The strain energy density along the principal fiber direction is typically higher than that perpendicular to the direction, reflecting the difference in the tissue's ability to store elastic energy in different directions. After calculating the strain energy density for each voxel, the results are stored as a three-dimensional scalar field. ,in Represents a spatial location vector.

[0111] For each feasible expansion path, a fast travel method is used to calculate the geodesic distance field. The path origin is set as the zero level set, and the distance values ​​of surrounding voxels are iteratively updated to ensure the distance field satisfies the equation of motion. The geodesic distance comprehensively considers Euclidean distance and the mechanical impedance of the structure the path traverses, resulting in a larger equivalent distance for paths passing through high-stiffness structures. Based on the calculated distance field, a threshold is set... Delineate the path influence domain This threshold is determined based on the expected effective range of the MSE device and is typically set to 1.5 to 2 times the path length. For voxels with a distance less than the threshold, their shortest distance to the path is calculated, and the influence weight is determined according to an exponential decay function. This gives higher weight to regions closer to the path.

[0112] When constructing an adaptive integral grid within the path influence domain, the grid density is dynamically adjusted based on the gradient magnitude of the direction-dependent strain energy density field. In regions where strain energy density changes drastically, the number of sampling points is increased by refining the grid to ensure integration accuracy. An octree data structure is used to manage the grid, starting with a coarse grid surrounding the entire influence domain and recursively subdividing cells to meet the refinement conditions. The refinement condition is based on the range of strain energy density within a cell; refinement is triggered when the range exceeds a certain percentage threshold of the global range. After multiple levels of refinement, an integral node is established at the center of the leaf node cell, and integration weights are assigned according to the cell size.

[0113] At each integration node, the spatial partial derivatives of the direction-dependent strain energy density field are calculated using a central difference scheme. Strain energy density values ​​are collected at adjacent positions before and after the node along the three coordinate axes, and the partial derivative components are approximated using difference quotients. , , The three partial derivative components are combined into a gradient vector. Spatial partial derivatives are calculated for each of the three components of the gradient vector, yielding nine second-order partial derivatives. The sum of the three diagonal terms constitutes the result of the divergence operator acting on the gradient, i.e., the Laplace operator acting on the strain energy density field. The result of this divergence operation characterizes the source and sink distribution of energy per unit volume and is defined as energy flux density. Positive values ​​indicate energy convergence, while negative values ​​indicate energy diffusion.

[0114] Extract the outer boundary surface of the path influence domain, which is composed of isosurfaces of the distance field. Triangularly mesh the boundary surface and calculate the unit outward normal vector at the centroid of each triangular element. The normal vector is obtained by normalizing the cross product of the coordinates of the three vertices of the facet element, and its direction points outward from the influence domain. At the centroid of each facet element, the energy flux density value is sampled, and the energy flux density is considered as the expression of the vector field at that point. The dot product of the energy flux density vector and the boundary normal vector is calculated to obtain the energy flux density passing through that facet element. The sign of this dot product reflects the energy flow direction; a positive value indicates energy flowing out of the influence domain, and a negative value indicates energy flowing in. The energy flux density of each facet element is multiplied by the area of ​​the facet element to obtain the total energy flux on that facet element. Then, the energy fluxes of all facet elements are summed to complete the surface integral along the boundary.

[0115] When performing volume integration within the path's influence domain, a previously constructed adaptive integration mesh is utilized. All leaf node cells in the mesh are traversed, and the direction-dependent strain energy density value is extracted at the integration node of each cell. This density is multiplied by the cell's volume and the influence weighting function value to obtain the node's contribution to the total volume integration. The contributions of all integration nodes are summed to obtain the total strain energy within the influence domain. This volume integration reflects the total elastic energy stored in the tissue surrounding the path due to the expansion of the MSE device. The energy flux obtained from the boundary surface integration is added to the internal strain energy obtained from the volume integration; according to the principle of energy conservation, the sum equals the total deformation energy contribution corresponding to the path. Considering the interaction between the path and the surrounding tissue, a coupling coefficient is introduced. The modified boundary flux term, which reflects the efficiency of energy transfer at the boundary, typically ranges from 0.6 to 0.9.

[0116] For multiple feasible expansion paths, the deformation energy associated with each path reflects the mechanical response strength of the organization when that path is chosen. To comprehensively evaluate the overall effect of all paths, a nonlinear aggregation function is employed. This function is designed to consider the interactions between paths, avoiding overestimation caused by simple linear superposition. An aggregation form based on Minkowski sums is used to aggregate the deformation energy of each path. Summing after exponentiation, then taking the first power. The root of the power. Parameters Controlling the degree of nonlinearity in aggregation, when When equal to 1, it degenerates into linear summation; when It approaches the maximum value as it approaches infinity. Based on experimental observations of the tissue's mechanical response, we set... A value between 2 and 3 indicates that paths with higher deformation energy dominate the aggregation result, while other paths still contribute to some extent. In the numerical calculation, the deformation energy value for each path... (in After performing exponentiation on the path index and summing the results, the total tissue deformation energy is finally obtained through radical calculation. This total amount serves as the objective function term in the subsequent optimization process. By minimizing the total deformation energy, it controls tissue damage and ensures that the expansion scheme of the MSE device maximizes the protection of surrounding physiological structures while achieving the corrective goal.

[0117] Figure 2 This is a flowchart illustrating the deformation energy calculation process of an MSE device based on a digital intelligent MSE automatic design method using biomechanical simulation, as described in an embodiment of the present invention.

[0118] In one alternative implementation,

[0119] The displacement field and strain field are obtained by iterative calculation using the spatial configuration parameters as boundary conditions. The cumulative deformation is obtained by time integration of the displacement field, and the iterative deformation energy is obtained by spatial integration of the strain field, including:

[0120] The position coordinates and expansion direction parameters of the MSE device are extracted from the spatial configuration parameters. Based on the position coordinates, the position of the point of action is calibrated in the three-dimensional morphological data and the force boundary conditions are constructed in combination with the expansion direction parameters. Fixed constraint regions are identified from the three-dimensional morphological data. Zero displacement constraints are applied to the fixed constraint regions to obtain displacement boundary conditions and combined with the force boundary conditions to obtain complete boundary conditions.

[0121] The three-dimensional morphological data is discretized using the finite element method to construct an overall stiffness matrix. The complete boundary conditions are applied to the overall stiffness matrix to obtain the displacement vectors of each node. Based on the displacement vectors, shape function interpolation is performed to obtain a continuous displacement field. The spatial partial derivatives of the displacement vectors of each voxel in the continuous displacement field are calculated to obtain the displacement gradient tensor. The displacement gradient tensor is symmetric to obtain the strain tensor and arranged to obtain the strain field. The displacement increment vectors of the continuous displacement field between adjacent time steps are calculated and accumulated along the time axis using the trapezoidal integral method to obtain the cumulative deformation.

[0122] Anisotropic elastic tensors are extracted from the mechanical response characteristics of the tissue and bilinearly contracted with the strain tensors of each voxel in the strain field to obtain the voxel strain energy density. The voxel strain energy density is then volume integraled in the spatial domain to obtain the iterative deformation energy.

[0123] After obtaining the spatial configuration parameters of the MSE device, the effectiveness of the configuration needs to be verified and the tissue response calculated through biomechanical simulation. The spatial configuration parameters include the three-dimensional position information and expansion direction information of the MSE device. These parameters need to be converted into boundary condition forms recognizable by finite element analysis. Specifically, the position coordinates are given as triplets in Cartesian coordinates, and the expansion direction parameter is represented as a unit direction vector. When locating the point of application of the MSE device in the three-dimensional morphological data, a nearest neighbor search algorithm is used to find the node with the smallest Euclidean distance to the given position coordinates in the mesh node set as the force loading point. The expansion direction parameter determines the direction of the applied force. Combined with the expansion force amplitude set based on clinical experience, a complete concentrated force boundary condition is constructed. This force boundary condition is expressed in the form of nodal load vectors, where the magnitude of the force is determined according to the design specifications of the MSE device, typically in the range of 50N to 200N.

[0124] Identifying fixed constraint regions relies on anatomical knowledge and clinical experience. In maxillary expansion scenarios, the skull base suture region, the zygomatic junction, and the bony structures at the base of the nasal cavity are typically considered relatively fixed constraint regions. These regions are identified in 3D morphological data through anatomical landmarks or obtained through bone mineral density thresholding. For all mesh nodes within the identified fixed constraint regions, three translational constraints with zero degrees of freedom are applied, meaning that nodes in that region are subject to these constraints. , , The displacement components in all three directions are forced to zero. The force boundary conditions and displacement boundary conditions together constitute the complete boundary conditions. The former describes the action of external loads, while the latter defines the motion boundary of the structure. Only by combining the two can the finite element equations have a unique solution.

[0125] Three-dimensional morphological data is discretized using tetrahedral mesh elements. The mesh density is adaptively adjusted based on geometric complexity, with localized refinement in areas of high curvature or high stress gradient. Each tetrahedral element contains four nodes, each with three translational degrees of freedom; therefore, the stiffness matrix of a single element is 12×12. The global stiffness matrix is ​​obtained by assembling the element stiffness matrices. The assembly process follows the correspondence between node numbers and degree-of-freedom numbers, accumulating the elements from each element's stiffness matrix into the corresponding positions in the global stiffness matrix. The global stiffness matrix has the following dimensions: ,in This represents the total number of grid nodes. This matrix exhibits sparsity and symmetric positive definiteness, making it suitable for solving using the conjugate gradient method or the preconditional conjugate gradient method.

[0126] The complete boundary conditions are applied to the global stiffness matrix, which is then modified to reflect the constraints. For displacement boundary conditions, a "set to 1" method is used: the diagonal elements of the row corresponding to the degree of freedom are set to 1, all other elements in that row are set to 0, and the right-hand side is set to the constraint displacement value (0 in this case). For force boundary conditions, the nodal load vectors are directly superimposed onto the corresponding positions of the global load vector. The modified linear equations are in the form of... ,in This is the corrected global stiffness matrix. Let be the nodal displacement vector to be determined. This is the corrected load vector. Solving this system of equations yields the displacement components of all nodes; these discrete nodal displacement values ​​constitute the fundamental data of the displacement field.

[0127] The construction of a continuous displacement field relies on shape function interpolation techniques. For tetrahedral elements, using linear shape functions, the displacement at any point within the element can be expressed as: ,in For the first The shape function of each node. Let be the displacement vector of this node. The coordinates are local coordinates within the element. Shape functions satisfy the property of taking a value of 1 at the corresponding node and a value of 0 at other nodes, and the sum of the four shape functions is always 1. Using this interpolation method, displacement values ​​can be calculated at any location within the element, thus obtaining a spatially continuous displacement field.

[0128] The strain field is calculated based on geometric equations under the assumption of small deformation. The spatial partial derivatives of the displacement field constitute the displacement gradient tensor, whose components are... ,in Indicates displacement at The directional components. The displacement gradient tensor contains information about both rigid body rotation and pure deformation, while the strain tensor only reflects pure deformation and is obtained through symmetry processing, i.e., the strain components. ,in, This represents the component of the displacement vector in the b-direction. For a tetrahedral linear element, the partial derivatives of the shape function with respect to spatial coordinates are constant within the element, therefore the strain field is uniformly distributed within the element. Arranging the strain tensors of all elements according to their spatial positions forms a discrete representation of the overall strain field.

[0129] The cumulative deformation is calculated using time-dimensional integration. The expansion process of the MSE device is gradual, typically lasting several weeks, with small expansion increments each day. In the numerical simulation, the entire expansion process is discretized into multiple time steps, each corresponding to one expansion increment. The displacement increment vector between adjacent time steps is defined as... superscript This indicates the time step number. The trapezoidal integration method accumulates the displacement increments, and the accumulated deformation is expressed as... ,in For the first The step size of each time step. Let be the displacement rate. Under the quasi-static loading assumption, the displacement rate can be approximated as the displacement increment divided by the time step. Therefore, the cumulative deformation is simplified to the weighted sum of the displacement increments at each time step. The cumulative deformation reflects the overall displacement response of the tissue throughout the expansion process and is an important indicator for evaluating the treatment effect.

[0130] Iterative deformation energy is calculated by combining the strain field and the material constitutive relation. The anisotropic elastic tensor stored in the micromechanical response characteristics is represented in the form of a fourth-order tensor, denoted as . Its weight Describing stress components With strain components The linear relationship between them. For anisotropic materials, the elastic tensor has 21 independent components, but under the assumption of orthogonal anisotropy, it can be simplified to 9 independent parameters. The voxel strain energy density is defined as... Using Einstein's summation convention, the summation is indicated by repeated indices. This expression is essentially a bilinear contraction operation of the strain tensor and the elasticity tensor, with the result being a scalar representing the elastic strain energy stored per unit volume.

[0131] The volume integral in the spatial domain sums the discrete voxel strain energy densities into the overall deformation energy. For a finite element mesh, the volume integral is transformed into a summation over all elements, and the iterative deformation energy is expressed as: ,in The total number of units, For the first Strain energy density of each element This represents the volume of the element. The volume of a tetrahedral element can be calculated using the determinant of the nodal coordinates. Iterative deformation energy reflects the overall mechanical response strength of the microstructure under the current configuration parameters; the smaller this value, the more uniform the microstructure deformation and the more reasonable the stress distribution.

[0132] In one alternative implementation,

[0133] Based on the difference between the cumulative deformation amount and the preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined, and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is then determined whether the iterative deformation energy satisfies the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions, including:

[0134] The deformation deviation is obtained by calculating the difference between the cumulative deformation and the preset physiological threshold. The deformation deviation is then judged by sign and marked as exceeding the limit. The difference between the iterative deformation energy and the iterative deformation energy corresponding to the previous iteration is calculated, and the rate of change of iterative deformation energy is calculated based on the difference.

[0135] Logarithmic transformation is performed on the absolute value of the deformation deviation to obtain the deviation scale, and reciprocal transformation is performed on the absolute value of the rate of change of the iterative deformation energy to obtain the convergence rate. Hyperbolic tangent mapping is performed on the deviation scale and the convergence rate to obtain the adaptive step size factor, and the adjustment step size is calculated by combining the preset benchmark step size.

[0136] The gradient vector is obtained by calculating the partial derivatives of the iterative deformation energy with respect to each parameter component in the spatial configuration parameters. The gradient vector is normalized by the second norm to determine the unit gradient direction. The correction direction coefficient is determined based on the over-limit state, and the correction increment vector is calculated by combining the unit gradient direction and the adjustment step size. The correction increment vector is superimposed on the spatial configuration parameters to obtain the correction configuration parameters. It is determined whether the absolute value of the rate of change of the iterative deformation energy is less than a preset convergence threshold. If it is less, the correction configuration parameters are output as design instructions.

[0137] After completing the iterative calculations of the displacement and strain fields, the spatial configuration parameters need to be adaptively adjusted to bring the design scheme closer to the optimal state. The cumulative deformation at the end of the current iteration cycle is extracted; this cumulative deformation reflects the overall displacement response of the target anatomical structure under the action of the MSE device. The cumulative deformation is compared with a preset physiological threshold, which is usually set based on clinical experience. For example, the safe range for unilateral expansion during maxillary suture expansion is 4 to 8 mm; exceeding this range may cause complications such as root resorption or bone fenestration. The difference between the cumulative deformation and the preset physiological threshold is calculated to obtain the deformation deviation. The deformation deviation can be positive or negative. A positive value indicates that the current cumulative deformation exceeds the physiological threshold, while a negative value indicates that the expected expansion target has not yet been reached. The sign of the deformation deviation is determined; when... The time marker indicates an over-limit state as "over-expansion". The time marker indicates an "insufficient expansion" state. The time stamp is marked as "compliant". The marking of the out-of-limit state will play a key role in determining the direction of subsequent corrections.

[0138] Evaluate the convergence characteristics of the iterative process. Record the iterative deformation energy of the current iteration cycle. And extract the iterative deformation energy corresponding to the previous iteration from historical data. Calculate the difference between the two. This difference reflects the changing trend of the tissue deformation energy across successive iterations. The rate of change of iterative deformation energy is calculated based on this difference. Rate of change The smaller the absolute value, the closer the iterative process is to a steady state, and the more gradual the energy change of the system. In practical applications, when Less than It is generally assumed that the iteration is close to convergence. It is important to note that the sign of the rate of change also has physical meaning: a negative value indicates that the system energy is decreasing, while a positive value indicates that the energy is increasing. The latter may indicate that the configuration parameters are adjusted in an inappropriate direction or that the step size is too large, leading to oscillations.

[0139] To achieve adaptive step size control, the deformation deviation and the rate of change of iterative deformation energy need to be converted into dimensionless adjustment factors. The absolute value of the deformation deviation... Perform a logarithmic transformation to obtain the deviation scale. ,in For reference deformation, 10% of a preset physiological threshold is typically used as the normalization benchmark. Logarithmic transformation can compress the influence of large deviations and prevent the step size from becoming uncontrollable due to over-response. The absolute value of the rate of change of the iterative deformation energy... Perform a reciprocal transformation to obtain the convergence rate. ,in This is a numerically stable term, with a value of approximately To prevent the denominator from being zero. Convergence rate quantity The larger the step size, the slower the iteration changes, and the more the step size should be increased to accelerate convergence; conversely, the smaller the step size, the more the step size should be decreased to ensure stability.

[0140] Deviation scale With convergence rate Perform hyperbolic tangent mapping to obtain the adaptive step size factor. ,in and These are the mapping weight coefficients, typically taken as... and This reflects a strategy that prioritizes bias correction over convergence acceleration. The hyperbolic tangent function maps the input to... The interval ensures that the step size factor is always positive and bounded. This is combined with a preset baseline step size. The adjustment step size is calculated. The baseline step size is set according to the problem size. For the maxillary suture expansion problem, the typical baseline step size is 0.05 mm, which corresponds to a device position or orientation adjustment of about 0.05 mm in each iteration.

[0141] After determining the adjustment step size, the correction direction needs to be calculated. Spatial configuration parameters typically include the installation location coordinates of the MSE device. , expansion constraint direction and optimal path length Multiple components. Calculate iterative deformation energy. The gradient vector is obtained by taking the partial derivatives with respect to each parameter component. The partial derivatives are calculated using the finite difference method, applying a small perturbation to each parameter component. The iterative deformation energy is recalculated by approximating the derivative using the forward difference or central difference formula. The gradient vector points in the direction of the fastest increase in iterative deformation energy, and its negative direction is the direction of the fastest decrease in energy.

[0142] Normalizing the gradient vector using the second norm yields the unit gradient direction. ,in This represents the Euclidean norm. Normalization eliminates the influence of gradient magnitude, making the correction direction determined solely by the geometric properties of energy changes. The correction direction coefficients are determined based on the overlimit state. When the over-limit state is "over-expansion", it is necessary to reduce the tissue deformation energy to reduce the intensity of mechanical stimulation. In this case, take... Correction is made along the negative gradient direction; when the over-limit state is "insufficient expansion", the tissue deformation energy needs to be increased to enhance the expansion effect, at which point... Correct along the positive gradient direction; when the out-of-limit state is the "acceptable state", take... Pause correction. Combine with unit gradient direction. and adjusting step size The corrected increment vector is calculated. .

[0143] The incremental vector will be corrected. Overlay to current space configuration parameters Obtain the corrected configuration parameters The revised configuration parameters need to be checked for physical constraints. For example, the installation position coordinates must not exceed the anatomical structure boundary, the expansion constraint direction must maintain the unit vector property, and the optimal path length must not be negative. If the revised parameters violate the constraints, a projection correction is performed: the position coordinates are projected into the feasible region, the direction vector is renormalized, and the path length is truncated to a non-negative range.

[0144] After parameter correction, determine whether the iteration process has converged. Calculate the absolute value of the rate of change of the iterative deformation energy. and with the preset convergence threshold Comparisons are made. The convergence threshold is typically set to... to The specific value is determined based on the required precision of the problem. If This indicates that the energy change between successive iterations is small enough that the system has reached quasi-static equilibrium, at which point the output configuration parameters are corrected. As design instructions, these instructions include parameters such as the precise installation location, expansion direction, and expected expansion amount of the MSE device, and can be directly used for clinical surgical planning or digital manufacturing. Then the configuration parameters will be corrected. As the initial configuration for the next iteration, the displacement field and strain field calculations are re-executed, and a new iteration cycle begins.

[0145] In practical applications, to prevent iterations from getting stuck in local optima or infinite loops, a maximum iteration limit is typically set, such as 50 to 100 iterations. If convergence is not achieved after reaching the maximum iteration limit, the current optimal configuration parameters are output and marked as "partially converged," prompting clinicians to conduct manual review. Furthermore, an iteration history monitoring mechanism can be introduced to record key indicators such as deformation deviation, iteration deformation energy, and step size factor for each iteration, and to plot convergence curves to aid in diagnosing the stability and effectiveness of the iteration process. Through the aforementioned adaptive correction strategy, the spatial configuration parameters of the MSE device can be efficiently optimized to the optimal state that meets physiological constraints and mechanical performance requirements while ensuring computational stability.

[0146] A second aspect of this invention provides a digital intelligent MSE automatic design system based on biomechanical simulation, comprising:

[0147] The mechanical analysis unit is used to acquire three-dimensional morphological data and physiological load data of the target anatomical structure, map the physiological load data to the mechanical bearing path corresponding to the three-dimensional morphological data to obtain a load distribution matrix, and perform tensor operations on the load distribution matrix and pre-stored tissue mechanical response characteristics to obtain a stress prediction field.

[0148] The parameter optimization unit is used to identify high stress concentration areas and low stiffness weak areas based on the stress prediction field, perform boundary tracking and extract the contour curve of the high stress concentration area, calculate the normal gradient along the contour curve and determine the expansion constraint direction, and minimize the stress distribution variance and the total structural deformation energy according to the expansion constraint direction and the low stiffness weak area, and determine the spatial configuration parameters of the MSE device.

[0149] The iteration involves a unit used to iteratively calculate the displacement field and strain field by using the spatial configuration parameters as boundary conditions. The displacement field is integrated over time to obtain the cumulative deformation, and the strain field is integrated spatially to obtain the iterative deformation energy. Based on the difference between the cumulative deformation and a preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is then determined whether the iterative deformation energy meets the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions.

[0150] A third aspect of the present invention provides an electronic device, comprising:

[0151] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0152] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0153] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0154] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A digital intelligent MSE automatic design method based on biomechanical simulation, characterized in that, include: The three-dimensional morphological data and physiological load data of the target anatomical structure are acquired. The physiological load data is mapped to the mechanical bearing path corresponding to the three-dimensional morphological data to obtain the load distribution matrix. Tensor operations are performed on the load distribution matrix and the pre-stored tissue mechanical response characteristics to obtain the stress prediction field. Based on the stress prediction field, high stress concentration areas and low stiffness weak areas are identified. Boundary tracking is performed on the high stress concentration areas and contour curves are extracted. The normal gradient is calculated along the contour curves and the expansion constraint direction is determined. Based on the expansion constraint direction and the low stiffness weak areas, the stress distribution variance and the total structural deformation energy are minimized and the spatial configuration parameters of the MSE device are determined. The spatial configuration parameters are used as boundary conditions to iteratively calculate the displacement field and strain field. The displacement field is integrated over time to obtain the cumulative deformation, and the strain field is integrated spatially to obtain the iterative deformation energy. Based on the difference between the cumulative deformation and a preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined, and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is determined whether the iterative deformation energy meets the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions.

2. The method according to claim 1, characterized in that, The process involves acquiring three-dimensional morphological data and physiological load data of the target anatomical structure, mapping the physiological load data onto the mechanical bearing path corresponding to the three-dimensional morphological data to obtain a load distribution matrix, and performing tensor operations on the load distribution matrix and pre-stored tissue mechanical response features to obtain a stress prediction field, including: The tomographic image sequence of the target anatomical structure is acquired by medical imaging scanning equipment. The tomographic image sequence is reconstructed in three dimensions to obtain three-dimensional morphological data and the bone region and soft tissue region are segmented and extracted. The stress data of the target anatomical structure during physiological activities is acquired to obtain physiological load data and the direction is decomposed to obtain normal load component and tangential load component. Feature points are extracted based on the boundary of the skeletal region and a skeletal topological connection relationship is constructed. The stiffness distribution field is calculated based on the thickness distribution of the soft tissue region. The mechanical bearing path is determined based on the skeletal topological connection relationship and the stiffness distribution field. The normal load component is projected along the normal direction of the mechanical bearing path and accumulated to obtain the normal load distribution. The tangential load component is projected along the tangential direction of the mechanical bearing path and accumulated to obtain the tangential load distribution. A load distribution matrix is ​​constructed based on the normal load distribution and the tangential load distribution. The elastic modulus and Poisson's ratio are extracted from the mechanical response characteristics of the tissue. A stiffness tensor is constructed based on the elastic modulus and Poisson's ratio. The load distribution matrix is ​​transposed and tensor contraction is performed with the stiffness tensor to obtain the initial stress tensor. The principal stress components and shear stress components of the initial stress tensor are calculated and a stress prediction field is constructed.

3. The method according to claim 1, characterized in that, Based on the stress prediction field, high stress concentration regions and low stiffness weak regions are identified. Boundary tracing and contour curve extraction are performed on the high stress concentration regions. The normal gradient is calculated along the contour curve, and the expansion constraint direction is determined, including: Morphological expansion and corrosion operations are performed on the stress prediction field to extract local extreme value regions. The stress mean of the local extreme value regions is calculated and compared with the global stress mean to determine high stress concentration regions. Stiffness feature values ​​of each voxel are extracted from the three-dimensional morphological data and a stiffness field is constructed. The regions in the stiffness field with stiffness feature values ​​lower than a preset stiffness threshold and the regions in the stress prediction field with stress values ​​lower than a preset stress threshold are spatially intersected to obtain low stiffness weak regions. Voxelization boundary extraction is performed on the high stress concentration region to obtain a boundary voxel set. A three-dimensional point cloud is constructed based on the boundary voxel set and surface fitting is performed to obtain a boundary surface. Isoparametric sampling is performed on the boundary surface to obtain a sampling point sequence and a contour curve is constructed along the connection direction. At each sampling point on the contour curve, the unit normal vector of the boundary surface is calculated. The directional derivative of the stress prediction field is calculated along the direction of the unit normal vector to obtain the normal stress gradient. The normal stress gradient is integrated along the contour curve, and the normal vector with the largest integral value is extracted as the principal expansion direction. The curvature tensor of the contour curve in three-dimensional space is calculated, and the principal direction corresponding to the minimum curvature is extracted as the secondary expansion direction. The expansion constraint direction is constructed based on the principal expansion direction and the secondary expansion direction.

4. The method according to claim 1, characterized in that, Based on the expansion constraint direction and the low-stiffness weak region, minimizing the stress distribution variance and the total structural deformation energy, and determining the spatial configuration parameters of the MSE device, includes: Extract the unit vector of the expansion constraint direction and construct the direction matrix. Determine the centroid coordinates of the low stiffness weak region and the high stress concentration region respectively and calculate the connection vector. Calculate the dot product of the connection vector and the direction matrix to obtain the direction correlation. Based on the direction correlation, determine the spatial position type of the low stiffness weak region relative to the expansion constraint direction. Based on the spatial location type, candidate configuration regions are delineated in the three-dimensional morphological data and a set of sampling points is generated. The shortest distance between each sampling point in the set of sampling points and the boundary of the high stress concentration region is calculated to obtain the distance field. Combined with the expansion constraint direction, a feasible expansion path is determined. The path integral of the stress prediction field is performed along the feasible expansion path to obtain the cumulative path stress value. Combined with the length of the feasible expansion path, the average path stress is determined. Based on the average path stress, the change in stress distribution variance is determined. The strain energy density is calculated by extracting the elastic modulus distribution of the three-dimensional morphological data. The strain energy density is integrated within the volume covered by the feasible expansion path to obtain the total tissue deformation energy. The comprehensive evaluation index is calculated by combining the change in the stress distribution variance and the directional correlation. The sampling point with the minimum comprehensive evaluation index is identified in the candidate configuration area, and its coordinates and direction vectors are extracted as spatial configuration parameters.

5. The method according to claim 4, characterized in that, Extracting the elastic modulus distribution from the three-dimensional morphological data to calculate the strain energy density, and integrating the strain energy density over the volume covered by the feasible expansion path to obtain the total tissue deformation energy, including: Multi-scale tissue layering structure is extracted from the three-dimensional morphological data and fiber orientation tensor reconstruction is performed on each layer. An anisotropic elastic tensor field is constructed based on the fiber orientation tensor and strain tensor components are calculated by combining the stress prediction field. Bilinear tensor contraction operation is performed on the strain tensor components and the anisotropic elastic tensor field to obtain the direction-dependent strain energy density field. The geodesic distance field of the feasible expansion path is calculated and the path influence domain is constructed. Within the path influence domain, the direction-dependent strain energy density field is non-uniformly weighted and an adaptive integral grid is constructed. At each grid node, the local gradient of the direction-dependent strain energy density field is calculated and divergence is performed to obtain the energy flux density. The boundary energy flux is obtained by performing an inner product operation between the energy flux density and the boundary normal vector of the path influence domain. The boundary energy flux is then subjected to surface integral along the boundary of the path influence domain and combined with the volume integral within the adaptive integral grid to obtain the path-related deformation energy. The path-related deformation energy corresponding to all feasible expansion paths is then nonlinearly aggregated to obtain the total tissue deformation energy.

6. The method according to claim 1, characterized in that, The displacement field and strain field are obtained by iterative calculation using the spatial configuration parameters as boundary conditions. The cumulative deformation is obtained by time integration of the displacement field, and the iterative deformation energy is obtained by spatial integration of the strain field, including: The position coordinates and expansion direction parameters of the MSE device are extracted from the spatial configuration parameters. Based on the position coordinates, the position of the point of action is calibrated in the three-dimensional morphological data and the force boundary conditions are constructed in combination with the expansion direction parameters. Fixed constraint regions are identified from the three-dimensional morphological data. Zero displacement constraints are applied to the fixed constraint regions to obtain displacement boundary conditions and combined with the force boundary conditions to obtain complete boundary conditions. The three-dimensional morphological data is discretized using the finite element method to construct an overall stiffness matrix. The complete boundary conditions are applied to the overall stiffness matrix to obtain the displacement vectors of each node. Based on the displacement vectors, shape function interpolation is performed to obtain a continuous displacement field. The spatial partial derivatives of the displacement vectors of each voxel in the continuous displacement field are calculated to obtain the displacement gradient tensor. The displacement gradient tensor is symmetric to obtain the strain tensor and arranged to obtain the strain field. The displacement increment vectors of the continuous displacement field between adjacent time steps are calculated and accumulated along the time axis using the trapezoidal integral method to obtain the cumulative deformation. Anisotropic elastic tensors are extracted from the mechanical response characteristics of the tissue and bilinearly contracted with the strain tensors of each voxel in the strain field to obtain the voxel strain energy density. The voxel strain energy density is then volume integraled in the spatial domain to obtain the iterative deformation energy.

7. The method according to claim 1, characterized in that, Based on the difference between the cumulative deformation amount and the preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined, and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is then determined whether the iterative deformation energy satisfies the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions, including: The deformation deviation is obtained by calculating the difference between the cumulative deformation and the preset physiological threshold. The deformation deviation is then judged by sign and marked as exceeding the limit. The difference between the iterative deformation energy and the iterative deformation energy corresponding to the previous iteration is calculated, and the rate of change of iterative deformation energy is calculated based on the difference. Logarithmic transformation is performed on the absolute value of the deformation deviation to obtain the deviation scale, and reciprocal transformation is performed on the absolute value of the rate of change of the iterative deformation energy to obtain the convergence rate. Hyperbolic tangent mapping is performed on the deviation scale and the convergence rate to obtain the adaptive step size factor, and the adjustment step size is calculated by combining the preset benchmark step size. The gradient vector is obtained by calculating the partial derivatives of the iterative deformation energy with respect to each parameter component in the spatial configuration parameters. The gradient vector is normalized by the second norm to determine the unit gradient direction. The correction direction coefficient is determined based on the over-limit state, and the correction increment vector is calculated by combining the unit gradient direction and the adjustment step size. The correction increment vector is superimposed on the spatial configuration parameters to obtain the correction configuration parameters. It is determined whether the absolute value of the rate of change of the iterative deformation energy is less than a preset convergence threshold. If it is less, the correction configuration parameters are output as design instructions.

8. A digital intelligent MSE automatic design system based on biomechanical simulation, used to implement the method of any one of claims 1-7, characterized in that, include: The mechanical analysis unit is used to acquire three-dimensional morphological data and physiological load data of the target anatomical structure, map the physiological load data to the mechanical bearing path corresponding to the three-dimensional morphological data to obtain a load distribution matrix, and perform tensor operations on the load distribution matrix and pre-stored tissue mechanical response characteristics to obtain a stress prediction field. The parameter optimization unit is used to identify high stress concentration areas and low stiffness weak areas based on the stress prediction field, perform boundary tracking and extract the contour curve of the high stress concentration area, calculate the normal gradient along the contour curve and determine the expansion constraint direction, and minimize the stress distribution variance and the total structural deformation energy according to the expansion constraint direction and the low stiffness weak area, and determine the spatial configuration parameters of the MSE device. The iteration involves a unit used to iteratively calculate the displacement field and strain field by using the spatial configuration parameters as boundary conditions. The displacement field is integrated over time to obtain the cumulative deformation, and the strain field is integrated spatially to obtain the iterative deformation energy. Based on the difference between the cumulative deformation and a preset physiological threshold, and combined with the rate of change of the iterative deformation energy between consecutive iterations, the adjustment step size is determined and the spatial configuration parameters are corrected along the gradient direction to obtain the corrected configuration parameters. It is then determined whether the iterative deformation energy meets the preset convergence condition. If it does, the corrected configuration parameters are output as design instructions.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images