Neurosurgical procedure planning three-dimensional modeling system with digital twin fusion

Abstract

This invention relates to the field of three-dimensional medical image processing and geometric modeling technology, specifically disclosing a three-dimensional modeling system for neurosurgical planning that integrates digital twins. The system acquires multimodal images of patients and inverts modulus distribution to construct a three-dimensional solid mesh. Using the modulus step interface as a criterion to mark traversed elements, the shape function is reconstructed as a superposition format of step and transition compensation, and a forced transition constraint is applied to obtain a mesh dataset with transition bandwidth attributes. Based on the deformation gradient, the stress oscillation index of each element is calculated, and a time-domain piecewise solution strategy is selected to eliminate spurious oscillations and obtain the initial displacement field. A global regularization solution is implemented using the sum of squared spatial modulus gradients as a constraint to obtain a smooth displacement field. The smooth displacement field is superimposed back onto the initial topology and iteratively fine-tuned by minimizing length deviation, outputting a biomechanically corrected three-dimensional geometric mesh. This invention realizes the complete construction from a static anatomical model to a digital twin with biomechanical evolution capabilities.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional medical image processing and geometric modeling technology, specifically to a three-dimensional modeling system for neurosurgical planning that integrates digital twins. Background Technology

[0002] With the rapid development of high-field magnetic resonance imaging, diffusion tensor imaging, and multimodal image fusion technologies, neurosurgical planning has gradually moved from traditional two-dimensional image interpretation to three-dimensional visualization and quantitative analysis. Currently, three-dimensional geometric modeling methods based on individualized medical imaging data have been applied to surgical planning for complex intracranial diseases such as brain tumors, vascular malformations, and epilepsy. In the standard three-dimensional modeling process, multimodal image data is first registered and fused. The boundary contours of lesions and key anatomical structures are extracted using threshold segmentation, region growing, or deep learning-based segmentation algorithms. Then, surface reconstruction algorithms such as the traveling cube are used to generate a surface mesh in the form of triangular facets. Finally, a tetrahedralization algorithm is used to convert the surface mesh into a tetrahedral solid mesh suitable for finite element analysis. This process can effectively reproduce the macroscopic morphology of various intracranial tissue regions at the geometric level, providing surgeons with an intuitive three-dimensional anatomical reference.

[0003] Existing 3D solid meshes suffer from a mismatch between topological partitioning and the continuous distribution of material parameters. This mismatch leads to abrupt changes in interlayer parameters, causing abnormal stress wave propagation and spurious oscillations. Consequently, the displacement field solution is distorted, and the mesh edge length is warped, making it impossible for digital twins to reliably reproduce the deformation and evolution process of brain tissue. This invention reconstructs the shape function system to generate transition bandwidths at abrupt changes to eliminate stress shielding. It then obtains a continuous displacement field through oscillation suppression and regularization smoothing. Finally, iterative fine-tuning of the mesh coordinates by minimizing length deviations achieves coordination and unity of the 3D geometric mesh in three dimensions: topological preservation, material continuity, and mechanical consistency. Summary of the Invention

[0004] The purpose of this invention is to provide a three-dimensional modeling system for neurosurgical planning that integrates digital twins, in order to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: A 3D modeling system for neurosurgical planning that integrates digital twins includes: The data acquisition and mesh construction module collects multimodal medical image data of patients and modulus distribution obtained from image grayscale inversion, and constructs a three-dimensional solid mesh carrying modulus layering information. The interface marking and shape function reconstruction module marks the crossed elements of the 3D solid mesh using the modulus step interface as the criterion. The shape function system of the reconstructed elements is a superposition format of step components and transition compensation components. Forced transition constraints are applied to the interlayer modulus jump to obtain a mesh dataset with transition bandwidth properties. The oscillation suppression and time-domain piecewise solution module calculates the stress oscillation index of each element based on the current deformation gradient for the mesh dataset, and selects the time-domain piecewise solution strategy based on the element stress oscillation index to obtain the initial displacement field that eliminates false oscillations. The global regularization smoothing module performs global regularization on the initial displacement field with the sum of squared spatial modulus gradients as a constraint to obtain a smooth displacement field under continuous modulus transition. The displacement field superposition and output module superimposes the smooth displacement field onto the initial topology of the three-dimensional solid mesh to obtain a biomechanically corrected three-dimensional geometric mesh.

[0006] As a further aspect of the present invention: the construction process of the three-dimensional solid mesh is as follows: Extract the extreme points of gray-level gradient at the boundary between gray matter and white matter in multimodal medical image data to form a three-dimensional seed point set; Using a three-dimensional seed point set as anchor points, the isosurface normal of the modulus distribution is projected onto the local tangent plane at various sub-points to generate an irregular sampling field covering the entire brain space; Using the modulus values ​​at each point in the sampling field as weights, the node coordinates of the three-dimensional solid mesh are driven to make fine adjustments along their normals, thereby obtaining a three-dimensional solid mesh with the modulus layer boundary aligned with the spatial position of the mesh nodes.

[0007] As a further aspect of the present invention: obtaining a grid dataset with transition bandwidth attributes specifically includes: Based on the unit volume fraction divided by the modulus step interface, the attribute region to which each unit belongs is quantitatively defined, and the units in the transition interval are marked to cross the boundary, generating a list of units to be reconstructed. For each unit in the list of reconstructed units, the modulus difference is used as a constraint condition. Directional degrees of freedom are configured on the enriched nodes so that the step transition only takes effect along the normal direction of the interface, and the configured shape function format is obtained. After obtaining the shape function format, a transition constraint is constructed along the tangential direction of the interface to restrict the free movement of each enriched node to the same direction as the transition constraint, thereby obtaining the configuration state of each enriched node after constraint. The cells corresponding to the configured enriched nodes are combined into a new dataset. The node coordinates and connection relationships of the original cells are merged to obtain a grid dataset with transition bandwidth attributes.

[0008] As a further aspect of the present invention: configuring directional degrees of freedom on the enriched nodes specifically includes: Determine the normal direction of the interface where each enriched node is located, and calibrate the additional displacement degrees of freedom at the enriched nodes along the normal direction; Based on the direction calibration results of the normal direction, the tangential component of the additional displacement degree of freedom is constrained to zero, and only the independent motion permission of the normal component is retained. The additional displacement degrees of freedom, after orientation calibration and component constraints, are embedded into the shape function construction of each enriched element to obtain a configured shape function format that only takes effect along the interface normal direction.

[0009] As a further aspect of the present invention: obtaining the initial displacement field to eliminate spurious oscillations specifically includes: Based on the change in deformation gradient of each element at the current moment and the deformation gradient at the previous moment, the stress oscillation index value of each element is calculated. The stress oscillation index values ​​of each element are compared with preset thresholds. Elements that exceed the preset thresholds are given explicit time advance permissions to dissipate their oscillation components, while elements that do not exceed the preset thresholds are given implicit equilibrium solution permissions to maintain their quasi-static paths. The nodal displacements of each element after explicit time progression and those after implicit equilibrium solution are merged according to their respective nodal affiliations to obtain the initial displacement field that eliminates spurious oscillations.

[0010] As a further aspect of the present invention: the calculation of the stress oscillation index value of each unit specifically includes: Obtain the deformation gradient of each element at the current moment and the deformation gradient at the previous moment, and decompose them into volume change tensor and shape change tensor respectively. Calculate the first magnitude of change of the volume change tensor in the time domain and the second magnitude of change of the shape change tensor in the time domain, respectively. The first and second variation amplitudes are weighted and combined, and the result of the weighted combination is used as the stress oscillation index value.

[0011] As a further aspect of the present invention: obtaining a smooth displacement field under continuous modulus transition specifically includes: At each node of the 3D solid mesh, the sum of squared spatial modulus gradients of the node is calculated based on the modulus values ​​of the neighboring elements, and used as the constraint weight values ​​of the node. Using the constraint weights of each node as coefficients, amplitude suppression is applied to the displacement components of the initial displacement field at each node to generate a locally suppressed displacement field. The initial displacement field is replaced by a locally suppressed displacement field, which is used as the displacement input for global regularization solution to obtain a smooth displacement field.

[0012] As a further aspect of the present invention: the generation of the locally suppressed displacement field specifically includes: Arrange the constraint weight values ​​of each node in order of magnitude, take the median value of the arrangement as the suppression benchmark, and determine the suppression amplitude of the displacement component of each node by the difference between the weight value of each node and the suppression benchmark. In the displacement components of each node in the initial displacement field, the corresponding displacement component value is subtracted from the suppression amplitude determined for each node to obtain the amplitude-suppressed displacement component value of each node. The displacement components after amplitude suppression at each node are recombined to form a locally suppressed displacement vector at the node. The displacement vectors of all nodes constitute a locally suppressed displacement field.

[0013] As a further aspect of the present invention: the biomechanically corrected three-dimensional geometric mesh is obtained as follows: Extract the initial coordinate values ​​of each node in the initial topology of the 3D solid mesh, and obtain the displacement vectors of each node in the smooth displacement field; The initial coordinate values ​​and displacement vectors are superimposed on each node to obtain the temporary coordinate values ​​for each node. Centered on the temporary coordinates of each node, the length of the line connecting each neighboring node to the node is compared with the original length of the line in the initial topology. The temporary coordinates are iteratively fine-tuned based on the comparison deviation to minimize the total length deviation of all lines. The coordinates after fine-tuning convergence are used as the correction coordinates. The initial coordinate values ​​of the corresponding nodes in the initial topology are replaced with the corrected coordinate values ​​of each node to obtain a biomechanically corrected three-dimensional geometric mesh.

[0014] The beneficial effects of this invention are: (1) This invention marks the tetrahedral elements that pass through the modulus step interface and reconstructs their shape function system, and applies a forced transition constraint at the interlayer modulus jump, so that the material parameter step jump that originally existed due to topological discretization is converted into a transition bandwidth with physical continuity. This invention makes the modulus gradient inside and between three-dimensional mesh elements no longer produce non-physical abrupt changes, avoids the reflection and abnormal transmission attenuation of stress waves caused by parameter step at the material interface, and improves the stress-strain response accuracy of the three-dimensional model when simulating large deformation of tissue.

[0015] (2) This invention introduces a stress oscillation index judgment strategy based on the time-domain change of deformation gradient. It assigns explicit dissipation or implicit equilibrium solution permissions according to the oscillation state of each element, eliminating false high-frequency oscillation components induced by the jump of element properties. At the same time, it applies global regularization smoothing to the initial displacement field with the sum of squares of spatial modulus gradient as a constraint, so that the final output displacement field meets the physical constraint of continuous material modulus transition in space. The node coordinates are fine-tuned by minimizing the length deviation iteration, ensuring that the relative deformation of the mesh edge length reaches the global minimum while retaining the initial topological connection relationship. This provides a three-dimensional geometric benchmark with continuous spatial coordinates, consistent topological relationship and smooth physical properties for subsequent intraoperative navigation and mechanical simulation. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a system block diagram of the present invention; Figure 2 This is a flowchart of the process for obtaining a grid dataset with transition bandwidth attributes in this invention; Figure 3 This is a flowchart of the process for obtaining the initial displacement field that eliminates spurious oscillations in this invention. Detailed Implementation

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

[0019] Please see Figure 1 As shown, this invention is a 3D modeling system for neurosurgical planning that integrates digital twins, comprising: The data acquisition and mesh construction module collects multimodal medical image data of patients and modulus distribution obtained from image grayscale inversion, and constructs a three-dimensional solid mesh carrying modulus layering information. The interface marking and shape function reconstruction module marks the crossed elements of the 3D solid mesh using the modulus step interface as the criterion. The shape function system of the reconstructed elements is a superposition format of step components and transition compensation components. Forced transition constraints are applied to the interlayer modulus jump to obtain a mesh dataset with transition bandwidth properties. The oscillation suppression and time-domain piecewise solution module calculates the stress oscillation index of each element based on the current deformation gradient for the mesh dataset, and selects the time-domain piecewise solution strategy based on the element stress oscillation index to obtain the initial displacement field that eliminates false oscillations. The global regularization smoothing module performs global regularization on the initial displacement field with the sum of squared spatial modulus gradients as a constraint to obtain a smooth displacement field under continuous modulus transition. The displacement field superposition and output module superimposes the smooth displacement field onto the initial topology of the three-dimensional solid mesh to obtain a biomechanically corrected three-dimensional geometric mesh.

[0020] In the data acquisition and mesh construction module, multimodal medical image data of patients and modulus distribution obtained from image grayscale inversion are acquired, and a three-dimensional solid mesh carrying modulus layering information is constructed, specifically including: In the process of constructing a three-dimensional solid mesh, the extreme points of gray-level gradient at the boundary between gray matter and white matter are extracted from multimodal medical image data.

[0021] Specifically, the multimodal medical imaging data includes high-resolution T1-weighted and T2-weighted images, acquired by magnetic resonance imaging equipment with a slice thickness of 0.5 mm.

[0022] For each voxel in the above image data, calculate the absolute value of the difference in grayscale value between it and the twenty-six spatially adjacent voxels, and take the maximum value of the absolute value as the grayscale gradient value of the voxel.

[0023] Sort all voxels from high to low according to their gray-level gradient values, and take the top 5% of voxels as candidate extrema. Then, remove isolated points from the candidate extrema. The removal rule is that if there are no other candidate extrema within a radius of 3 voxels centered on a candidate extrema, then the point is removed. The remaining candidate extrema form a 3D seed point set.

[0024] Using a three-dimensional seed point set as anchor points, the isosurface normal of the modulus distribution is projected onto the local tangent plane at various sub-points to generate an irregular sampling field covering the entire brain space.

[0025] Modulus distribution is a three-dimensional scalar field composed of the elastic modulus values ​​of each voxel obtained by inverting the gray values ​​of an image.

[0026] For each seed point, firstly, the gradient direction of the modulus scalar field at that seed point is calculated; this gradient direction is the normal to the isosurface. Simultaneously, a plane perpendicular to the isosurface normal is constructed through that seed point; this is called the local tangent plane of that seed point.

[0027] The projection direction is obtained by projecting the isosurface normal onto the local tangent plane.

[0028] Along the projection direction, three voxel positions are taken on each side of the seed point as sampling points, so that the sampling points generated by each seed point cover all spatial positions in the whole brain tissue region, and all sampling points constitute an irregular sampling field.

[0029] Using the modulus values ​​at each point in the sampling field as weights, the node coordinates of the 3D solid mesh are finely adjusted along their normals to obtain a 3D solid mesh with the modulus layer boundary aligned with the spatial position of the mesh nodes. Specifically, for each node in the 3D solid mesh, all sampling points are searched within a radius of 1.5 voxels centered on that node. The modulus values ​​at each of the found sampling points are summed and averaged. This average value is then used as the representative modulus value for that node.

[0030] The normal of a node is taken as the average direction of the node's position at the initial moment pointing to all sampling points in its neighborhood.

[0031] Along the normal direction, the quotient obtained by dividing the modulus representative value by the maximum modulus value of the whole brain region is used as the step size coefficient. The step size coefficient is multiplied by the length of the voxel diagonal as the movement step size, so that the node moves along the normal direction by a movement step size. After completing a position fine adjustment for all nodes, the above fine adjustment operation is repeated three times to obtain a three-dimensional solid mesh with the modulus layer boundary aligned with the spatial position of the mesh nodes.

[0032] Please see Figure 2 As shown, in the interface marking and shape function reconstruction module, for a 3D solid mesh, the elements that are traversed are marked using the modulus step interface as the criterion. The shape function system of the reconstructed elements is a superposition format of step components and transition compensation components. A forced transition constraint is applied to the interlayer modulus jump to obtain a mesh dataset with transition bandwidth properties, specifically including: In the process of obtaining a grid dataset with transition bandwidth attributes, the attribute region to which each cell belongs is quantitatively defined based on the cell volume fraction segmented by the modulus step interface, and boundary crossing marks are applied to cells within the transition interval to generate a list of cells to be reconstructed.

[0033] The modulus step interface is the set of spatial locations in a three-dimensional solid mesh where the modulus values ​​of adjacent elements change abruptly.

[0034] For each tetrahedral element in a 3D solid mesh, the fraction of the element's total volume that is located on one side of the interface after the element is cut by the modulus step interface is calculated and is called the volume fraction.

[0035] If the volume fraction of a cell is between 0.3 and 0.7, the cell is defined as being within the transition range and marked as a boundary-crossing cell. If the volume fraction is less than 0.3, the cell belongs to a homogeneous region on one side of the interface. If the volume fraction is greater than 0.7, the cell belongs to a homogeneous region on the other side of the interface. After traversing all cells, all cells marked as boundary-crossing cells are summarized to generate a list of cells to be reconstructed.

[0036] For each element in the list of reconstructed elements, the modulus difference is used as a constraint condition. directional degrees of freedom are configured on the enriched nodes so that the step transition only takes effect along the normal direction of the interface, thus obtaining the configured shape function format.

[0037] Among them, enriched nodes refer to the nodes located on the modulus step interface among all the nodes contained in each unit in the list of units to be reconstructed.

[0038] The modulus difference refers to the difference in average modulus values ​​between the regions on both sides of the boundary cell associated with the node.

[0039] For each enriched node, an additional displacement degree of freedom is added at that node, independent of the node's original three spatial displacement degrees of freedom. The activation status of this additional degree of freedom is determined by the magnitude of the modulus difference; it is activated when the modulus difference is greater than 0.1 MPa, and deactivated otherwise.

[0040] For enriched nodes with additional degrees of freedom enabled, configure the additional degree of freedom to contribute displacement only along the interface normal direction.

[0041] Determine the normal direction of the interface where each enriched node is located, and calibrate the direction of the additional displacement degrees of freedom at the enriched node along the normal direction. Specifically, for each enriched node, take one element on each side of the interface from the boundary elements associated with the node, calculate the modulus gradient direction of these two elements respectively, and take the average direction of the two gradient directions as the normal direction of the interface where the enriched node is located.

[0042] The direction of the additional displacement degree of freedom at the enriched node is calibrated to be consistent with the normal direction.

[0043] Based on the orientation calibration results of the normal direction, the tangential component of the additional displacement degree of freedom is constrained to zero, and only the independent motion permission of the normal component is retained. That is, a local coordinate system is established at the enriched node, with one coordinate axis coinciding with the normal direction, and the other two coordinate axes located within the tangential plane of the interface. The component values ​​of the additional displacement degree of freedom in the two coordinate axis directions within the tangential plane are fixed to zero, and only the component value in the normal coordinate axis direction is retained as a variable.

[0044] The additional displacement degrees of freedom, after orientation calibration and component constraints, are embedded into the shape function construction of each enriched element. This involves superimposing the original standard shape function of each enriched element with the shape function contribution corresponding to the additional displacement degrees of freedom, so that the contribution only acts in the cross-interface normal direction, thus obtaining the configured shape function format.

[0045] After obtaining the above shape function format, a transition constraint is constructed along the tangential direction of the interface to restrict the free movement of each enriched node to the same direction as the transition constraint, thereby obtaining the configuration state of each enriched node after constraint.

[0046] The tangential direction of the interface is any direction within the plane perpendicular to the aforementioned normal direction. A constraint is applied in this tangential direction, stipulating that the rate of change of the additional displacement degrees of freedom of each enriched node in the tangential direction is zero. The elements corresponding to the configured enriched nodes are grouped into a new dataset, and the node coordinates and connectivity relationships of the original elements are merged to obtain a mesh dataset integrating the transition bandwidth attribute. Specifically: the original node coordinate values ​​of all tetrahedral elements and the connectivity relationships between nodes are retained unchanged, and the configured shape function format is written into the dataset as a new attribute of each enriched element, forming a mesh dataset integrating the transition bandwidth attribute.

[0047] Please see Figure 3 As shown, in the oscillation suppression and time-domain piecewise solution module, for the mesh dataset, the stress oscillation index of each element is calculated based on the current deformation gradient, and a time-domain piecewise solution strategy is selected based on the element stress oscillation index to obtain the initial displacement field that eliminates spurious oscillations. Specifically, this includes: In the process of obtaining the initial displacement field to eliminate spurious oscillations, the stress oscillation index value of each element is calculated based on the change amplitude of the deformation gradient of each element at the current moment and the deformation gradient at the previous moment.

[0048] The deformation gradient is a second-order physical quantity that describes the degree of deformation of an element from its initial configuration to its current configuration. It is determined by the relationship between the partial derivatives of the current coordinates and the initial coordinates of each node in the 3D solid mesh.

[0049] For each tetrahedral element, at the end of each time increment step, the deformation gradient value of the element at the current time is calculated based on the difference between the current coordinates of the four nodes of the element and the coordinates of the previous time step. Retrieve the deformation gradient value of the cell at the end of the previous time increment step from the storage cell, and use it as the deformation gradient of the previous time step.

[0050] After obtaining the deformation gradient at the current moment and the deformation gradient at the previous moment, decompose them into volume change tensor and shape change tensor, respectively.

[0051] The volume change tensor is characterized by the determinant of the deformation gradient, which is equal to the ratio of the current element volume to the initial element volume. The shape change tensor is characterized by the skew portion remaining after subtracting the volume change portion from the deformation gradient. Its magnitude is measured by the equivalent shear strain, which is defined as the square root of the sum of the squares of the components of the deviatoric strain tensor.

[0052] Calculate the first magnitude of the volume change tensor in the time domain and the second magnitude of the shape change tensor in the time domain, respectively.

[0053] The first change range is calculated using the following mathematical formula: ; in, This represents the determinant value of the deformation gradient at the current moment. This represents the determinant value of the deformation gradient at the previous moment. Let be the determinant of the deformation gradient at the initial moment, and take a value of 1. In this formula, when Calculations are performed if the value is not zero; if it is zero, the cell is skipped and marked as an exception. First change magnitude. It characterizes the relative fluctuation of volume change between adjacent time points; the larger the value, the more intense the volume oscillation.

[0054] The second change range is calculated using the following mathematical formula: ; in, The equivalent shear strain value of the shape change tensor at the current moment. This represents the equivalent shear strain value of the shape change tensor at the previous moment. The reference strain value is set to 0.01. and These represent the maximum and minimum equivalent shear strains recorded throughout the simulation process, initially set to 0.1 and 0 respectively, and updated after each calculation. When the value is zero, this item takes the value of 1 to avoid division by zero. Second variation range. It characterizes the intensity of shape changes between adjacent time points.

[0055] After obtaining the first and second variation amplitudes, the two are weighted and combined, and the result of the weighted combination is used as the stress oscillation index value.

[0056] The weighting coefficients of the weighted combination are: the first change amplitude is weighted at 0.6, the second change amplitude is weighted at 0.4, the two are multiplied and added together to obtain the stress oscillation index value of the unit at the current moment.

[0057] The stress oscillation index values ​​of each unit are compared with the preset threshold.

[0058] The preset threshold value is 0.05, which is determined based on statistical observations of the stable and unstable boundaries in multiple brain tissue simulation cases.

[0059] When the stress oscillation index of a certain element is greater than 0.05, it is determined that the element has significant spurious oscillation components. The element is then granted explicit time-progression permission, that is, displacement updates are performed using the central difference explicit scheme within the time increment step to dissipate its oscillation components. When the stress oscillation index is less than or equal to 0.05, the element is determined to be in a quasi-static equilibrium state. The element is then granted implicit equilibrium solution permission, that is, displacements are solved using the Newton iteration implicit scheme within the time increment step to maintain its quasi-static path.

[0060] The node displacements of each element after explicit time advancement and those after implicit equilibrium solution are merged according to their respective element affiliations. Specifically, for each node, all elements to which it belongs are collected, and the displacement contributions obtained by the node from the explicitly advanced elements and from the implicitly solved elements are recorded. The two are then weighted and averaged according to the number of elements to obtain the composite displacement value of the node.

[0061] After the synthesis of all nodes is completed, the synthesized displacement values ​​of each node are combined into a complete initial displacement field, which is the preliminary displacement result after eliminating spurious oscillations.

[0062] In the global regularization smoothing module, the initial displacement field is subjected to global regularization calculation with the sum of squared gradients of the spatial modulus as a constraint to obtain a smooth displacement field under continuous modulus transition, specifically including: In obtaining a smooth displacement field under continuous modulus transition, at each node of the 3D solid mesh, the sum of squares of the spatial modulus gradient of that node is calculated based on the modulus values ​​of its neighboring elements, and this value is used as the constraint weight value of that node. For each node in the 3D solid mesh, all tetrahedral elements to which that node belongs are determined, and this set is called the neighboring element set of that node.

[0063] For each cell in the neighborhood cell set, obtain the modulus values ​​at the four vertices of the cell, and calculate the modulus change rate of the cell in the three orthogonal directions in three-dimensional space. That is, subtract the modulus values ​​of the four vertices of the cell from each other and divide by the spatial distance between the corresponding vertices. Take the sum of the squares of the change rates in the three directions as the squared value of the modulus gradient of the cell.

[0064] The sum of the squared values ​​of the modulus gradients of all cells in the neighborhood cell set is the sum of the squared values ​​of the spatial modulus gradients at that node.

[0065] This value is used as the constraint weight value for the node. The larger the constraint weight value, the more drastic the spatial change in the modulus of the region where the node is located.

[0066] Using the constraint weights of each node as coefficients, amplitude suppression is applied to the displacement components of the initial displacement field at each node to generate a locally suppressed displacement field.

[0067] Each node in the initial displacement field contains displacement components in three directions, corresponding to the three orthogonal coordinate axes in three-dimensional space.

[0068] For each node, the constraint weight value of that node is taken as the suppression coefficient. The displacement component values ​​of that node in the three directions are multiplied by the ratio of the suppression coefficient to the maximum value of the constraint weight values ​​of all nodes. This ratio is a decimal between 0 and 1. This ratio is used as the retention ratio of the displacement component of that node.

[0069] Subtracting the retention ratio from 1 yields the suppression ratio. Multiplying the original displacement component value by the suppression ratio and then subtracting it from the original value completes the amplitude suppression of the displacement component at that node.

[0070] After performing the above operations on all nodes, each node obtains the values ​​of the three displacement components after amplitude suppression.

[0071] The process of generating the locally suppressed displacement field described above can be further refined as follows: the constraint weight values ​​of each node are arranged from low to high in terms of numerical value, and the value in the middle of the arrangement is taken as the suppression benchmark.

[0072] If the total number of nodes is even, the average of the two values ​​in the middle position of the permutation is taken as the suppression benchmark.

[0073] For each node, calculate the difference between the constraint weight value of the node and the suppression benchmark. If the difference is greater than zero, the ratio obtained by dividing the difference by the maximum difference between the constraint weight values ​​of all nodes and the suppression benchmark is used as the suppression amplitude of the displacement component of the node. If the difference is less than or equal to zero, the suppression amplitude of the node is set to zero.

[0074] In the displacement components of each node in the initial displacement field, according to the suppression amplitude determined for each node, the original displacement component values ​​of the node in the three directions are subtracted from the product of the original value and the suppression amplitude respectively to obtain the amplitude-suppressed displacement component values ​​of the node in the three directions.

[0075] The amplitude-suppressed displacement component values ​​of each node in the three directions are recombined to form a locally suppressed displacement vector at that node, which contains component values ​​in the three directions.

[0076] The displacement vectors of all nodes together constitute a locally suppressed displacement field.

[0077] The initial displacement field is replaced by the locally suppressed displacement field, which is used as the displacement input for global regularization solution to obtain the smooth displacement field.

[0078] The specific operation of global regularization is as follows: the displacement vectors of each node in the locally suppressed displacement field are used as the initial guess values, the sum of the squared values ​​of the spatial modulus gradients of all nodes is used as the regularization constraint term, and the constraint term is superimposed on the total energy functional of the displacement field.

[0079] During the solution process, the coefficient of the constraint term is set to 0.1. This coefficient is predetermined based on the principle of balancing the smoothness of the displacement field with the degree of data fitting.

[0080] The total energy functional is minimized by iterative solution. In each iteration, the constraint term value is calculated based on the current displacement field and compared with the previous iteration step. The iteration is terminated when the change in the constraint term value between two consecutive iterations is less than 0.001. The displacement field at the time of stopping the iteration is the smooth displacement field.

[0081] In the displacement field superposition and output module, the smoothed displacement field is superimposed onto the initial topology of the 3D solid mesh to obtain a biomechanically corrected 3D geometric mesh, specifically including: In the process of obtaining the biomechanically corrected three-dimensional geometric mesh, the initial coordinate values ​​of each node in the initial topology of the three-dimensional solid mesh are extracted, and the displacement vectors of each node in the smooth displacement field are obtained.

[0082] In the initial topology of the 3D solid mesh, each node stores the initial coordinate value of the node in 3D space, which includes the specific values ​​in the three coordinate axes.

[0083] In the smooth displacement field, a displacement vector is stored for each node, which contains the displacement component values ​​of that node in the three coordinate axis directions.

[0084] Extract the initial coordinate values ​​and displacement vectors of each node according to its corresponding node number.

[0085] The initial coordinate values ​​and displacement vectors are superimposed on each node to obtain temporary coordinate values ​​for each node.

[0086] For each node, the initial coordinate values ​​of the node in each of the three coordinate axes are added to the displacement components of the displacement vector in the corresponding three coordinate axes. The resulting new values ​​in the three coordinate axes are the temporary coordinate values ​​of the node.

[0087] After performing the above vector superposition on all nodes, all nodes obtain their own temporary coordinate values.

[0088] After obtaining the temporary coordinate values ​​of each node, the length of the line connecting each neighboring node to the node is compared with the original length of the line in the initial topology, centered on the temporary coordinate values ​​of each node. The temporary coordinate values ​​are iteratively fine-tuned according to the comparison deviation to minimize the total length deviation of all lines, and the coordinate values ​​after fine-tuning convergence are used as the correction coordinate values.

[0089] The initial topology of the 3D solid mesh records the original length value of each line consisting of two nodes.

[0090] For each node, determine all its neighboring nodes that are directly connected to it through grid edges in the initial topology, forming the set of neighboring nodes for that node.

[0091] Using the temporary coordinates of the node as the center, calculate the spatial distance between the node and each neighboring node in its set of neighboring nodes to obtain the current connection length.

[0092] Subtract the current connection length from the original length of the corresponding connection in the initial topology to obtain the length deviation of the connection.

[0093] If the length deviation is positive, it indicates that the current connection is stretched; if it is negative, it indicates that the current connection is compressed.

[0094] The specific operation of iterative fine-tuning is as follows: For each node, collect the length deviation of the connection between it and all neighboring nodes, calculate the arithmetic mean of these deviations, and use half of this average value as the coordinate adjustment amount of the node in the current iteration step.

[0095] If the average deviation is positive, the temporary coordinate value of the node is moved by the distance of the coordinate adjustment in the direction pointing to the center position of the node and all its neighboring nodes; if the average deviation is negative, the coordinate adjustment is moved by the distance of the coordinate adjustment in the opposite direction.

[0096] After simultaneously adjusting the coordinates of all nodes, recalculate the total length deviation of all connected lines.

[0097] The method for calculating the total length deviation is as follows: add the absolute value of the difference between the current length and the original length of each line to obtain the total value.

[0098] Repeat the above iterative fine-tuning until the change in the sum of the length deviations obtained from two consecutive iterations is less than one ten-thousandth of the original sum of the lengths of all the connections. At this point, the iteration is considered to have converged, the fine-tuning is stopped, and the coordinate values ​​of each node at the time of convergence are used as the correction coordinate values.

[0099] The initial coordinate values ​​of the corresponding nodes in the initial topology are replaced with the corrected coordinate values ​​of each node to obtain a biomechanically corrected three-dimensional geometric mesh.

[0100] That is, while preserving the connection relationship of all nodes and the composition of elements in the initial topology of the three-dimensional solid mesh, the spatial position coordinates of each node are updated from the initial coordinate values ​​to the corrected coordinate values ​​to obtain a three-dimensional geometric mesh with a spatial configuration after biomechanical deformation.

[0101] The working principle of this invention is as follows: Multimodal medical image data of patients and the modulus distribution obtained from inversion are acquired to construct a three-dimensional solid mesh carrying modulus hierarchical information. The crossed elements are marked using the modulus step interface as a criterion. The shape function system of the reconstructed elements is a superposition format of step components and transition compensation components, and a forced transition constraint is applied to obtain a mesh dataset with transition bandwidth attributes. The stress oscillation index of each element is calculated based on the current deformation gradient. An explicit or implicit time-domain piecewise solution strategy is selected to eliminate spurious oscillations and obtain the initial displacement field. The initial displacement field is then subjected to global regularization using the sum of squared spatial modulus gradients as a constraint to obtain a smooth displacement field under continuous modulus transition. Finally, the smooth displacement field is superimposed back onto the initial topology of the three-dimensional solid mesh and fine-tuned iteratively by minimizing the length deviation to output a biomechanically corrected three-dimensional geometric mesh, thus realizing the complete construction from a static anatomical model to a digital twin with biomechanical evolution capabilities.

[0102] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A three-dimensional modeling system for neurosurgical planning integrating digital twins, characterized in that, include: The data acquisition and mesh construction module collects multimodal medical image data of patients and modulus distribution obtained from image grayscale inversion, and constructs a three-dimensional solid mesh carrying modulus layering information. The interface marking and shape function reconstruction module marks the crossed elements of the 3D solid mesh using the modulus step interface as the criterion. The shape function system of the reconstructed elements is a superposition format of step components and transition compensation components. Forced transition constraints are applied to the interlayer modulus jump to obtain a mesh dataset with transition bandwidth properties. The oscillation suppression and time-domain piecewise solution module calculates the stress oscillation index of each element based on the current deformation gradient for the mesh dataset, and selects the time-domain piecewise solution strategy based on the element stress oscillation index to obtain the initial displacement field that eliminates false oscillations. The global regularization smoothing module performs global regularization on the initial displacement field with the sum of squared spatial modulus gradients as a constraint to obtain a smooth displacement field under continuous modulus transition. The displacement field superposition and output module superimposes the smooth displacement field onto the initial topology of the three-dimensional solid mesh to obtain a biomechanically corrected three-dimensional geometric mesh.

2. The three-dimensional modeling system for neurosurgical planning integrating digital twins as described in claim 1, characterized in that, The process of constructing the three-dimensional solid mesh is as follows: Extract the extreme points of gray-level gradient at the boundary between gray matter and white matter in multimodal medical image data to form a three-dimensional seed point set; Using a three-dimensional seed point set as anchor points, the isosurface normal of the modulus distribution is projected onto the local tangent plane at various sub-points to generate an irregular sampling field covering the entire brain space; Using the modulus values ​​at each point in the sampling field as weights, the node coordinates of the three-dimensional solid mesh are driven to make fine adjustments along their normals, thereby obtaining a three-dimensional solid mesh with the modulus layer boundary aligned with the spatial position of the mesh nodes.

3. The three-dimensional modeling system for neurosurgical planning integrating digital twins as described in claim 1, characterized in that, Obtaining the grid dataset with transition bandwidth attributes specifically includes: Based on the unit volume fraction divided by the modulus step interface, the attribute region to which each unit belongs is quantitatively defined, and the units in the transition interval are marked to cross the boundary, generating a list of units to be reconstructed. For each unit in the list of reconstructed units, the modulus difference is used as a constraint condition. Directional degrees of freedom are configured on the enriched nodes so that the step transition only takes effect along the normal direction of the interface, and the configured shape function format is obtained. After obtaining the shape function format, a transition constraint is constructed along the tangential direction of the interface to restrict the free movement of each enriched node to the same direction as the transition constraint, thereby obtaining the configuration state of each enriched node after constraint. The cells corresponding to the configured enriched nodes are combined into a new dataset. The node coordinates and connection relationships of the original cells are merged to obtain a grid dataset with transition bandwidth attributes.

4. The three-dimensional modeling system for neurosurgical planning integrating digital twins according to claim 3, characterized in that, The configuration of directional degrees of freedom on enriched nodes specifically includes: Determine the normal direction of the interface where each enriched node is located, and calibrate the additional displacement degrees of freedom at the enriched nodes along the normal direction; Based on the direction calibration results of the normal direction, the tangential component of the additional displacement degree of freedom is constrained to zero, and only the independent motion permission of the normal component is retained. The additional displacement degrees of freedom, after orientation calibration and component constraints, are embedded into the shape function construction of each enriched element to obtain a configured shape function format that only takes effect along the interface normal direction.

5. The three-dimensional modeling system for neurosurgical planning integrating digital twins according to claim 1, characterized in that, Obtaining the initial displacement field to eliminate spurious oscillations specifically includes: Based on the change in deformation gradient of each element at the current moment and the deformation gradient at the previous moment, the stress oscillation index value of each element is calculated. The stress oscillation index values ​​of each element are compared with preset thresholds. Elements that exceed the preset thresholds are given explicit time advance permissions to dissipate their oscillation components, while elements that do not exceed the preset thresholds are given implicit equilibrium solution permissions to maintain their quasi-static paths. The nodal displacements of each element after explicit time progression and those after implicit equilibrium solution are merged according to their respective nodal affiliations to obtain the initial displacement field that eliminates spurious oscillations.

6. The three-dimensional modeling system for neurosurgical planning integrating digital twins according to claim 5, characterized in that, The calculation of the stress oscillation index value of each unit specifically includes: Obtain the deformation gradient of each element at the current moment and the deformation gradient at the previous moment, and decompose them into volume change tensor and shape change tensor respectively. Calculate the first magnitude of change of the volume change tensor in the time domain and the second magnitude of change of the shape change tensor in the time domain, respectively. The first and second variation amplitudes are weighted and combined, and the result of the weighted combination is used as the stress oscillation index value.

7. The three-dimensional modeling system for neurosurgical planning integrating digital twins according to claim 1, characterized in that, The process of obtaining a smooth displacement field under continuous modulus transition specifically includes: At each node of the 3D solid mesh, the sum of squared spatial modulus gradients of the node is calculated based on the modulus values ​​of the neighboring elements, and used as the constraint weight values ​​of the node. Using the constraint weights of each node as coefficients, amplitude suppression is applied to the displacement components of the initial displacement field at each node to generate a locally suppressed displacement field. The initial displacement field is replaced by a locally suppressed displacement field, which is used as the displacement input for global regularization solution to obtain a smooth displacement field.

8. The three-dimensional modeling system for neurosurgical planning integrating digital twins according to claim 7, characterized in that, The generation of the locally suppressed displacement field specifically includes: Arrange the constraint weight values ​​of each node in order of magnitude, take the median value of the arrangement as the suppression benchmark, and determine the suppression amplitude of the displacement component of each node by the difference between the weight value of each node and the suppression benchmark. In the displacement components of each node in the initial displacement field, the corresponding displacement component value is subtracted from the suppression amplitude determined for each node to obtain the amplitude-suppressed displacement component value of each node. The displacement components after amplitude suppression at each node are recombined to form a locally suppressed displacement vector at the node. The displacement vectors of all nodes constitute a locally suppressed displacement field.

9. The three-dimensional modeling system for neurosurgical planning integrating digital twins according to claim 1, characterized in that, The biomechanically corrected three-dimensional geometric mesh is obtained as follows: Extract the initial coordinate values ​​of each node in the initial topology of the 3D solid mesh, and obtain the displacement vectors of each node in the smooth displacement field; The initial coordinate values ​​and displacement vectors are superimposed on each node to obtain the temporary coordinate values ​​for each node. Centered on the temporary coordinates of each node, the length of the connection between each neighboring node and the node is compared with the original length of the connection in the initial topology. The temporary coordinates are iteratively fine-tuned based on the comparison deviation to minimize the total length deviation of all connections. The coordinates after fine-tuning convergence are used as the correction coordinates. The initial coordinates of the corresponding nodes in the initial topology are replaced with the correction coordinates of each node to obtain the biomechanically corrected three-dimensional geometric mesh.

Images