A method and system for segmenting and reconstructing three-dimensional medical image data

Abstract

The application relates to the technical field of medical image analysis, and discloses a three-dimensional medical image data segmentation and reconstruction method and system. The application effectively solves the contradiction between medical image denoising and edge preservation by using an anisotropic diffusion technology guided by a structure tensor field, significantly enhances the signal-to-noise ratio of a weak boundary, automatically removes vascular and skeletal interference by using eigenvalue analysis of a Hessian matrix, accurately locates a lesion topological core, eliminates subjective errors of manual interaction, gives the level set evolution an automatic truncation capability at a weak boundary adhesion place by using a curvature reverse brake mechanism, completely prevents segmentation leakage by using physical surface tension, and finally realizes mechanical balance of a grid vertex at sub-voxel accuracy by extremum solving of a variational energy functional, generates a three-dimensional entity model with geometric high precision and topological smoothness, and meets strict requirements of precision medicine on surgical planning and biomechanical analysis.

Description

Technical Field

[0001] This application relates to the field of medical image analysis technology, and in particular to a method and system for segmentation and reconstruction of three-dimensional medical image data. Background Technology

[0002] With the rapid development of computer-aided diagnosis and treatment technologies, three-dimensional reconstruction of medical images has become a crucial element in the planning of precision surgery and radiotherapy. In clinical practice, especially in the diagnosis and treatment of abdominal and pelvic diseases, medical professionals heavily rely on high-resolution computed tomography (CT) data to analyze the anisotropic spatial relationships between lesions and surrounding tissues. Unlike high-contrast tissues such as bones, soft tissues in the abdominal and pelvic cavities, such as the liver, pancreas, intestines, and intricate vascular networks, possess extremely complex physiological structures and very high density overlap. High-precision three-dimensional models reconstructed from two-dimensional image sequences can intuitively present the geometry, size, and spatial adjacency of lesions, providing surgeons with quantifiable anatomical information to assist in developing optimal preoperative surgical pathways and precise intraoperative navigation, thereby significantly reducing surgical trauma and improving patient prognosis. Therefore, the accuracy of image segmentation directly determines the realism of subsequent three-dimensional model reconstruction and the safety of clinical decisions, serving as a core foundation for realizing digital healthcare.

[0003] However, generating high-quality 3D models of soft tissue lesions faces significant technical challenges in actual data processing. The core difficulty lies in the inherent partial volume effect of computed tomography (CT) imaging and the diverse density distribution of soft tissues. In contrast-enhanced CT images of the abdomen or pelvis, the CT values ​​of malignant tumors are often highly similar to those of adjacent muscle tissue, blood vessel walls, or interstitial edema zones. The grayscale differences sometimes fluctuate within a very narrow range, and the two are often physically closely adhered. Existing image segmentation techniques, such as region growing algorithms or global threshold-based segmentation methods, primarily rely on the similarity of pixel grayscale values ​​and spatial connectivity for region extraction. When encountering situations where lesions are physically adhered to surrounding normal organs or have blurred boundaries, these algorithms lack effective boundary recognition mechanisms based on texture features or geometric topology. The algorithm evolution process is prone to penetrating weak boundaries with indistinct grayscale changes, incorrectly identifying non-target areas as lesion tissue. This segmentation overflow and oversegmentation phenomenon directly results in the generated segmentation mask volume being significantly larger than the actual anatomical structure, causing severe 3D morphological distortion. To correct such segmentation errors, clinical practice currently relies primarily on manual trimming and contour drawing layer by layer. This is not only inefficient and time-consuming, but also makes it difficult to guarantee the repeatability and objectivity of the results. If such 3D models with overflow errors are applied to surgical guide design or radiotherapy target delineation, it is highly likely to induce doctors to misjudge the resection range, damaging normal functional tissues, or causing deviations in radiation dose distribution, resulting in irreversible iatrogenic harm to patients. Therefore, how to achieve adaptive truncation of weak boundaries in complex scenarios with low contrast and tissue adhesion, ensuring the geometric and topological accuracy of segmentation and reconstruction, has become a pressing technical challenge in this field. Summary of the Invention

[0004] This application proposes a method and system for segmentation and reconstruction of three-dimensional medical image data to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, this application adopts the following technical solution: a method for segmentation and reconstruction of three-dimensional medical image data, comprising the following steps:

[0006] Step S1: Receive three-dimensional medical volumetric image data, calculate the gradient information of voxels in the three-dimensional medical volumetric image data to construct a structure tensor field, perform eigenvalue decomposition on the structure tensor field to obtain structure tensor eigenvalues, construct an anisotropic diffusion impedance field based on the structure tensor eigenvalues, substitute the anisotropic diffusion impedance field into the diffusion equation to perform directional smoothing processing on the three-dimensional medical volumetric image data, and output a preprocessed image volume containing an enhanced gradient vector field.

[0007] Step S2: Receive the preprocessed image volume output in step S1, construct a Hessian matrix for the voxels in the preprocessed image volume, calculate the eigenvalues ​​of the Hessian matrix, remove tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix, and extract regions exhibiting gray-scale energy convergence characteristics as the initial seed point set.

[0008] Step S3: Receive the initial seed point set output in step S2, construct an initial level set function with the initial seed point set as the zero level set region, construct a curvature inverse braking term, construct a gradient driving term based on the enhanced gradient vector field in step S1, establish a level set evolution equation composed of the curvature inverse braking term and the gradient driving term, perform iterative evolution on the initial level set function using the level set evolution equation, stop local evolution when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function makes the curvature inverse braking term greater than the gradient driving term, and output the evolved level set function.

[0009] Specifically, the gradient driving term is decoupled into two independent parts: a vector advection traction term and a scalar expansion thrust term. The inner product of the enhanced gradient vector field and the gradient vector of the initial level set function is calculated, and the inner product is defined as the vector advection traction term. The product of the preset expansion pressure coefficient and the edge indicator scalar function is calculated, and the product is defined as the scalar expansion thrust term.

[0010] An adaptive surface tension mechanism is constructed based on Yang-Laplace's law. The normalized divergence of the gradient vector of the initial level set function is calculated. The normalized divergence of the gradient vector is defined as the local mean curvature of the zero level set interface. A surface tension coefficient is introduced and set as a function of the edge indicator scalar function. This allows the surface tension coefficient to automatically increase at weak boundaries where the edge indicator scalar function has a large value. The product of the surface tension coefficient, the edge indicator scalar function, and the local mean curvature is calculated and defined as the curvature inverse braking term.

[0011] A partial differential evolution equation following the Hamiltonian Jacobian form is established, and the evolution rate of the partial differential evolution equation is determined by the algebraic sum of the vector advection traction term, the scalar expansion thrust term, and the curvature reverse braking term.

[0012] Step S4: Receive the evolved level set function output from step S3. Extract the initial mesh vertices and initial vertex connections from the zero level set interface of the evolved level set function. Calculate the data attachment energy term for the initial mesh vertices based on the enhanced gradient vector field from step S1. Calculate the elastic deformation energy term for the initial mesh vertices based on the initial vertex connections. Construct a variational energy functional consisting of a weighted sum of the data attachment energy term and the elastic deformation energy term. Calculate the optimized spatial coordinates of the initial mesh vertices by solving for the minimum value of the variational energy functional. Generate a 3D mesh model based on the optimized spatial coordinates.

[0013] Specifically, for each initial grid vertex, the interpolation gradient vector of the initial grid vertex at continuous physical space coordinates is calculated using the trilinear interpolation operator; the negative value of the square of the magnitude of the interpolation gradient vector is defined as the data attachment energy term;

[0014] The mean curvature vector of the initial mesh vertices is calculated using the discrete Laplace Beltrami operator. The mean curvature vector is orthogonally decomposed into a tangential relaxation component and a normal smoothing component. Tangential tension is constructed based on the tangential relaxation component, and normal pressure is constructed based on the normal smoothing component. A local adaptive stiffness coefficient is constructed, which is set to be inversely proportional to the local modulus of the enhanced gradient vector field. The local adaptive stiffness coefficient is set to approach zero at strong edges with large local modulus of the enhanced gradient vector field to close the normal pressure, and is set to increase at weak edges with small local modulus of the enhanced gradient vector field to maintain the shape using surface tension. The algebraic sum of the tangential tension and the normal pressure weighted by the local adaptive stiffness coefficient is calculated, and the algebraic sum is defined as the elastic deformation energy term.

[0015] Furthermore, in step S1, the specific operation of calculating the gradient information of voxels in the three-dimensional medical volumetric image data to construct the structure tensor field is as follows:

[0016] A physical observation scale is defined, which is used to filter out quantum noise in 3D medical volumetric image data. A Gaussian kernel function that conforms to the physical observation scale is used to perform 3D convolution smoothing on the 3D medical volumetric image data to generate a continuous differentiable physical field.

[0017] Calculate the gradient vector of each voxel in a continuously differentiable physical field in three-dimensional space. The gradient vector indicates the direction of gray-scale change of the continuously differentiable physical field at the physical observation scale.

[0018] Calculate the outer product matrix of the gradient vector and its transpose. The outer product matrix describes the tensor product of the gradient directions at a single point.

[0019] Define a physical integral scale, which determines the range of the physical support domain of local geometric perception. Use a Gaussian weighting function that conforms to the physical integral scale to perform a weighted average calculation on the outer product matrix within the range of the physical support domain to generate a second-order symmetric positive semi-definite matrix describing the statistical distribution of local gradient energy. The second-order symmetric positive semi-definite matrix is ​​the structure tensor field.

[0020] The specific operation of eigenvalue decomposition to obtain structural tensor eigenvalues ​​is as follows: the structural tensor field is decomposed into three non-negative structural tensor eigenvalues ​​and three mutually orthogonal structural tensor eigenvectors. The structural tensor eigenvalue with the largest value among the three structural tensor eigenvalues ​​represents the structural strength along the normal direction of the anatomical structure, and the two structural tensor eigenvalues ​​with the smaller values ​​among the three structural tensor eigenvalues ​​represent the structural strength along the tangential direction of the anatomical structure. By comparing the numerical differences among the three structural tensor eigenvalues, flat regions, planar structural regions, and tubular structural regions in three-dimensional medical volumetric imaging data are distinguished.

[0021] Furthermore, in step S1, the specific operation of constructing the anisotropic diffusion impedance field based on the structural tensor eigenvalues ​​is as follows:

[0022] An anisotropic diffusion tensor is constructed, which is the anisotropic diffusion impedance field. The eigenvectors of the anisotropic diffusion tensor are set to be consistent with the eigenvectors of the structure tensor of the structure tensor field. The eigenvalues ​​of the anisotropic diffusion tensor are defined as heat flux conductivity.

[0023] Set a tissue contrast sensitivity threshold, which is used to distinguish random thermal noise from solid anatomical boundaries at the energy level.

[0024] Calculate the square of the difference between the largest and smallest structure tensor eigenvalues ​​among the three structure tensor eigenvalues, and define the square of the difference as the coherence index.

[0025] For the heat flux conductivity along the anatomical structure normal, a negative exponential decay function based on the coherence index and the tissue contrast sensitivity threshold is used for calculation. This makes the heat flux conductivity along the anatomical structure normal approach zero when the coherence index is greater than the tissue contrast sensitivity threshold to form an adiabatic barrier, and makes the heat flux conductivity along the anatomical structure normal approach one when the coherence index is less than the tissue contrast sensitivity threshold to allow heat flow.

[0026] For the heat flux conductivity along the tangential direction of the anatomical structure, the value of the heat flux conductivity along the tangential direction of the anatomical structure is kept constant at one to allow omnidirectional diffusion;

[0027] The specific operation of performing directional smoothing on 3D medical volumetric image data by substituting the anisotropic diffusion impedance field into the diffusion equation is as follows: establish a diffusion evolution equation that follows the principle of non-equilibrium thermodynamics, use the anisotropic diffusion impedance field to correct the diffusion flux in the diffusion evolution equation, and as physical time progresses, use natural diffusion to eliminate noise in the flat areas of the 3D medical volumetric image data, and use normal blocking to maintain the boundary in the edge areas of the 3D medical volumetric image data, until the diffusion evolution equation reaches a quasi-steady state, and output a preprocessed image volume containing an enhanced gradient vector field.

[0028] Furthermore, in step S2, the specific operation of constructing a Hessian matrix for the voxels in the preprocessed image volume and calculating the eigenvalues ​​of the Hessian matrix is ​​as follows:

[0029] A discrete physical observation aperture sequence containing multiple Gaussian smoothing scale parameters with different values ​​is defined, and a corresponding three-dimensional Gaussian kernel function is constructed for each Gaussian smoothing scale parameter in the discrete physical observation aperture sequence.

[0030] A three-dimensional convolution operation is performed on the volume of the preprocessed image using a three-dimensional Gaussian kernel function to generate a multi-scale spatial response field. The square of the Gaussian smoothing scale parameter is calculated, and the square of the Gaussian smoothing scale parameter is defined as the scale normalization factor. The multi-scale spatial response field is multiplied and weighted using the scale normalization factor to generate a scale-normalized multi-scale spatial response field. The second-order partial derivative of each voxel in the scale-normalized multi-scale spatial response field in the Cartesian coordinate system is calculated.

[0031] A 3x3 second-order symmetric matrix describing the curvature of the local energy surface is constructed using second-order partial derivatives. The 3x3 second-order symmetric matrix is ​​defined as the Hessian matrix. Eigenvalue decomposition is performed on the Hessian matrix to obtain three Hessian matrix eigenvalues ​​characterizing the principal curvature. The three Hessian matrix eigenvalues ​​are then sorted in ascending order according to their absolute values.

[0032] The Hessian matrix eigenvalue with the smallest absolute value among the three Hessian matrix eigenvalues ​​is defined as the first principal curvature along the structural axis, the Hessian matrix eigenvalue with the largest absolute value among the three Hessian matrix eigenvalues ​​is defined as the third principal curvature along the cross-sectional direction, and the Hessian matrix eigenvalue with the middle absolute value among the three Hessian matrix eigenvalues ​​is defined as the second principal curvature.

[0033] Furthermore, in step S2, the specific operation of removing tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix and extracting regions exhibiting gray-level energy convergence characteristics as the initial seed point set is as follows:

[0034] Calculate the absolute value of the first principal curvature and define the absolute value of the first principal curvature as the minimum curvature intensity term. The minimum curvature intensity term is used to ensure that the voxel has significant curvature in all directions to eliminate planar structures.

[0035] Calculate the ratio of the absolute value of the first principal curvature to the absolute value of the third principal curvature. Define the ratio of the absolute value of the first principal curvature to the absolute value of the third principal curvature as the geometric isotropic factor. The geometric isotropic factor is used to quantify the degree to which the local geometric structure approximates a sphere in order to eliminate tubular structures.

[0036] Define a polarity gate operator, which is used to detect the sign consistency of the three Hessian matrix eigenvalues. When it is necessary to extract bright lesions, configure the polarity gate operator to output a non-zero value when all three Hessian matrix eigenvalues ​​are negative. When it is necessary to extract dark lesions, configure the polarity gate operator to output a non-zero value when all three Hessian matrix eigenvalues ​​are positive.

[0037] The product of the minimum curvature intensity term, the geometric isotropic factor, and the polarity gate operator is calculated, and the product is defined as the structural saliency response value of the voxel under the Gaussian smoothing scale parameter.

[0038] The maximum value of the structural significance response of the voxel under all Gaussian smoothing scale parameters is selected, and the maximum value is defined as the multi-scale convergence probability field. The multi-scale convergence probability field represents the probability that the voxel belongs to the core of the lesion.

[0039] Calculate the non-zero mean and non-zero standard deviation of the multi-scale convergent probability field, calculate the sum of the products of the non-zero mean, non-zero standard deviation and confidence factor, and define the sum as the adaptive truncation threshold.

[0040] Regions with values ​​greater than the adaptive truncation threshold in the multi-scale convergence probability field are extracted, and a binarized mask is generated. The binarized mask is defined as the initial seed point set.

[0041] Furthermore, in step S3, the initial level set function is constructed with the initial seed point set as the zero level set region, and the curvature inverse braking term is constructed. The specific operation of constructing the gradient driving term based on the enhanced gradient vector field in step S1 is as follows:

[0042] Calculate the shortest Euclidean distance from each voxel in 3D space to the edge of the initial seed point set. Construct a symbolic distance field based on the shortest Euclidean distance. The symbolic distance field is defined to take negative values ​​inside the initial seed point set, positive values ​​outside the initial seed point set, and zero values ​​at the edge of the initial seed point set. Define the symbolic distance field as the initial level set function.

[0043] Define a narrow band region where the initial level set function is close to zero. Update the initial level set function only within the narrow band region. Calculate the vector magnitude of the enhanced gradient vector field output in step S1. Construct an edge indicator scalar function based on the vector magnitude of the enhanced gradient vector field. Set the edge indicator scalar function to approach zero at strong edges where the vector magnitude of the enhanced gradient vector field is large, and set the edge indicator scalar function to approach one at flat regions where the vector magnitude of the enhanced gradient vector field is small.

[0044] Among them, the vector advection traction term drives the zero-level set interface to adhere to the anatomical structure boundary along the direction of the enhanced gradient vector field, while the scalar expansion thrust term simulates the fluid injection pressure driving the zero-level set interface to expand outward.

[0045] Furthermore, in step S3, an evolution equation for the level set is established, consisting of a curvature inverse braking term and a gradient driving term. This equation is used to iteratively evolve the initial level set function. Local evolution stops when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function causes the curvature inverse braking term to be greater than the gradient driving term. The specific operation for outputting the evolved level set function is as follows:

[0046] A partial differential evolution equation following the Hamiltonian Jacobi form is established. The partial differential evolution equation describes the dynamic deformation process of the initial level set function with physical time. A regularized Dirac function is introduced into the partial differential evolution equation, and the evolution force is restricted to the physical location of the zero level set interface by the regularized Dirac function.

[0047] In the iterative evolution process, for the flat internal region of the anatomical structure, the scalar expansion thrust term dominates, driving the zero-level set interface to grow outward.

[0048] For the strong edge regions of the anatomical structure, the edge indicator scalar function value approaches zero, causing the evolution rate to approach zero, locking the zero-level set interface on the boundary of the anatomical structure. For the weak boundary adhesion regions of the anatomical structure, the narrow channel causes the local average curvature value to increase sharply, so that the contraction force generated by the curvature reverse braking term exceeds the expansion force generated by the scalar expansion thrust term, causing the local evolution rate to become negative, forcing the zero-level set interface to automatically shrink back and stop growing.

[0049] The spatial derivative of the initial level set function is discretized using the central difference method, and the temporal derivative of the initial level set function is discretized using the forward Euler method. The time step is dynamically adjusted according to the numerical stability condition, which requires that the time step be less than the minimum value of the ratio of the spatial grid size to the maximum evolution rate.

[0050] The iterative evolution is continuously executed until the partial differential evolution equation reaches a steady state or the total energy change rate is lower than a preset threshold. The final steady-state scalar field is output as the evolved level set function, and the zero level set interface of the evolved level set function is extracted as the surface of the segmented three-dimensional anatomical structure.

[0051] Furthermore, in step S4, the evolved level set function output from step S3 is received, and the specific operation of extracting the initial mesh vertices and initial vertex connection relationships from the zero level set interface of the evolved level set function is as follows:

[0052] The zero level set interface in the evolved level set function output by step S3 is parsed using the moving cube isosurface extraction algorithm. The zero level set interface is parsed to obtain the discrete initial mesh vertex set and the initial vertex connection relationship. The initial mesh vertex set contains the initial mesh vertex coordinate information defined in three-dimensional space, and the initial vertex connection relationship contains the triangle mesh list that defines the physical topological adjacency attribute between the initial mesh vertices.

[0053] The specific operation for calculating the data attachment energy term for the initial grid vertices based on step S1 of enhancing the gradient vector field is as follows: construct a trilinear interpolation operator for the enhanced gradient vector field defined on the discrete voxel grid points output by step S1, and use the trilinear interpolation operator to establish a mapping relationship from the discrete data space to the continuous physical space, so as to transform the discrete enhanced gradient vector field into a continuous gradient potential energy density field.

[0054] The data-attached energy term physically constructs a continuous gravitational potential well. The energy valley of the gravitational potential well strictly corresponds to the gradient maximum position of the real physical boundary of the anatomical structure. The gravitational potential well generates a gradient adsorption force pointing towards the gradient maximum position. The gradient adsorption force drives the initial mesh vertices to break through the voxel resolution limit. The gradient adsorption force drives the initial mesh vertices to automatically slide towards the position of the local gradient energy maximum in the continuous physical space with sub-voxel precision.

[0055] Furthermore, in step S4, the elastic deformation energy term for the initial mesh vertices is calculated based on the initial vertex connection relationship. A variational energy functional is constructed, consisting of a weighted sum of the data attachment energy term and the elastic deformation energy term. The optimized spatial coordinates of the initial mesh vertices are calculated by solving for the minimum value of the variational energy functional. The specific operation of generating a 3D mesh model based on the optimized spatial coordinates is as follows:

[0056] The average curvature vector of the initial mesh vertices is calculated using the discrete Laplace Beltrami operator, and an area-weighted normal vector is introduced to prevent numerical singularities when the initial mesh vertices degenerate.

[0057] Tangential tension optimizes the shape quality and vertex distribution uniformity of the triangular mesh in the initial vertex connection relationship, while strictly maintaining the volume of the 3D mesh model;

[0058] Normal pressure eliminates high-frequency noise on the initial mesh vertex surface;

[0059] A variational energy functional containing data attachment energy and elastic deformation energy is constructed. The variational energy functional is iteratively solved using the semi-implicit Euler method until the maximum displacement of the initial grid vertex is less than the preset thermodynamic convergence threshold.

[0060] The output is the optimized spatial coordinates of the initial mesh vertex when it reaches mechanical equilibrium. The initial vertex connection relationship remains unchanged. The coordinates of the initial mesh vertex are updated to the optimized spatial coordinates to generate a 3D mesh model with sub-voxel level geometric accuracy and smooth topology.

[0061] A segmentation and reconstruction system for three-dimensional medical image data includes: an image preprocessing module, a seed point extraction module, a level set evolution module, and a mesh model generation module, wherein;

[0062] The image preprocessing module receives three-dimensional medical volumetric image data, calculates the gradient information of voxels in the three-dimensional medical volumetric image data to construct a structure tensor field, performs eigenvalue decomposition on the structure tensor field to obtain structure tensor eigenvalues, constructs an anisotropic diffusion impedance field based on the structure tensor eigenvalues, substitutes the anisotropic diffusion impedance field into the diffusion equation to perform directional smoothing processing on the three-dimensional medical volumetric image data, and outputs a preprocessed image volume containing an enhanced gradient vector field.

[0063] The seed point extraction module receives the preprocessed image volume output by the image preprocessing module, constructs a Hessian matrix for the voxels in the preprocessed image volume, calculates the eigenvalues ​​of the Hessian matrix, removes tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix, and extracts regions exhibiting gray-scale energy convergence characteristics as the initial seed point set.

[0064] The level set evolution module receives the initial seed point set output by the seed point extraction module, constructs an initial level set function with the initial seed point set as the zero level set region, constructs a curvature inverse braking term, constructs a gradient driving term based on the gradient vector field enhanced by the image preprocessing module, establishes a level set evolution equation composed of the curvature inverse braking term and the gradient driving term, performs iterative evolution on the initial level set function using the level set evolution equation, stops local evolution when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function makes the curvature inverse braking term greater than the gradient driving term, and outputs the evolved level set function.

[0065] The mesh model generation module receives the evolved level set function output by the level set evolution module, extracts the initial mesh vertices and initial vertex connections from the zero level set interface of the evolved level set function, calculates the data attachment energy term for the initial mesh vertices based on the enhanced gradient vector field of the image preprocessing module, calculates the elastic deformation energy term for the initial mesh vertices based on the initial vertex connections, constructs a variational energy functional composed of the weighted sum of the data attachment energy term and the elastic deformation energy term, calculates the optimized spatial coordinates of the initial mesh vertices by solving the minimum value of the variational energy functional, and generates a three-dimensional mesh model based on the optimized spatial coordinates.

[0066] The beneficial effects of this invention are as follows:

[0067] This invention effectively resolves the contradiction between denoising and edge preservation in medical images through anisotropic diffusion technology guided by structural tensor fields, significantly enhancing the signal-to-noise ratio of weak boundaries. Utilizing Hessian matrix eigenvalue analysis, it automatically eliminates vascular and skeletal interference, accurately locates the topological core of lesions, and eliminates subjective errors from manual interaction. A curvature-reverse braking mechanism endows the level set evolution with automatic truncation capabilities at weak boundary adhesion points. Physical surface tension is used to completely prevent segmentation leakage. Finally, through extreme value solving of variational energy functionals, mechanical equilibrium of mesh vertices at sub-voxel precision is achieved, generating a three-dimensional solid model with both high geometric accuracy and topological compliance, meeting the stringent requirements of precision medicine for surgical planning and biomechanical analysis. Attached Figure Description

[0068] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort:

[0069] Figure 1 This is a flowchart of the method of the present invention;

[0070] Figure 2 This is a system framework diagram of the present invention. Detailed Implementation

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

[0072] Example 1

[0073] like Figure 1 As shown, this invention discloses a method for segmenting and reconstructing three-dimensional medical image data, comprising the following steps:

[0074] Step S1: Receive three-dimensional medical volumetric image data, calculate the gradient information of voxels in the three-dimensional medical volumetric image data to construct a structure tensor field, perform eigenvalue decomposition on the structure tensor field to obtain structure tensor eigenvalues, construct an anisotropic diffusion impedance field based on the structure tensor eigenvalues, substitute the anisotropic diffusion impedance field into the diffusion equation to perform directional smoothing processing on the three-dimensional medical volumetric image data, and output a preprocessed image volume containing an enhanced gradient vector field.

[0075] Specifically, in step S1, the specific operation of calculating the gradient information of voxels in the three-dimensional medical volumetric image data to construct the structure tensor field is as follows:

[0076] A physical observation scale is defined, which is used to filter out quantum noise in 3D medical volumetric image data. A Gaussian kernel function that conforms to the physical observation scale is used to perform 3D convolution smoothing on the 3D medical volumetric image data to generate a continuous differentiable physical field.

[0077] In this embodiment, the three-dimensional medical volumetric image data comes from a computed tomography (CT) scanner. The physical unit of the voxel values ​​in the three-dimensional medical volumetric image data is the Henle unit, with a range of -1000 Henle units to 3000 Henle units.

[0078] Considering the Poisson distribution noise caused by photon statistical fluctuations during computed tomography imaging, directly performing differential operations on three-dimensional medical volumetric image data would lead to amplification of high-frequency noise. Therefore, a physical observation scale is introduced.

[0079] In this embodiment, the specific value of the physical observation scale is selected as 0.7 mm. The basis for choosing 0.7 mm as the physical observation scale is that the inter-slice resolution of conventional abdominal computed tomography scans is usually 0.625 mm to 1.25 mm. If the physical observation scale is less than 0.5 mm, quantum noise cannot be effectively suppressed; if the physical observation scale is greater than 1.0 mm, it will cause the edges of small lesions to be blurred. This step uses three-dimensional convolution operation to convolve the three-dimensional medical volumetric image data with a Gaussian kernel function whose standard deviation is equal to the physical observation scale, and outputs a continuously differentiable physical field. The continuously differentiable physical field establishes the differentiable manifold structure of the three-dimensional medical volumetric image data in Euclidean space.

[0080] Calculate the gradient vector of each voxel in a continuously differentiable physical field in three-dimensional space. The gradient vector indicates the direction of gray-scale change of the continuously differentiable physical field at the physical observation scale.

[0081] This step calculates the partial derivatives of each voxel in the continuously differentiable physical field along the three orthogonal axes of the Cartesian coordinate system by performing a central difference operation on the continuously differentiable physical field. These partial derivatives are then combined to form a three-dimensional gradient vector. The magnitude of the gradient vector characterizes the local contrast intensity at the physical observation scale, and the direction of the gradient vector points to the direction of the normal that increases the fastest in the Heinrich's unit.

[0082] Calculate the outer product matrix of the gradient vector and its transpose. The outer product matrix describes the tensor product of the gradient directions at a single point.

[0083] This step generates a 3x3 symmetric matrix by performing matrix multiplication on the gradient vector of each voxel, calculating the product of the gradient vector column vector and the gradient vector transpose row vector.

[0084] Define a physical integral scale, which determines the physical support domain of local geometric perception. Use a Gaussian weighting function that conforms to the physical integral scale to perform a weighted average calculation on the outer product matrix within the physical support domain to generate a second-order symmetric positive semi-definite matrix describing the statistical distribution of local gradient energy. The second-order symmetric positive semi-definite matrix is ​​the structure tensor field.

[0085] In this embodiment, the specific value of the physical integration scale is selected as 1.5 mm. The basis for selecting 1.5 mm as the physical integration scale is that the average radius of curvature of blood vessels and tumor edges in the human abdomen is usually greater than 1.5 mm. Setting the physical integration scale to 1.5 mm can ensure that the structural tensor field captures the geometric orientation of the macroscopic anatomical structure during the integration process, while suppressing local gradient fluctuations caused by tissue texture.

[0086] This step uses the current voxel as the center and a Gaussian weighted function with a standard deviation equal to the physical integral scale to perform a weighted average summation of all outer product matrices within the physical support domain, outputting a structure tensor field. Each element in the structure tensor field represents the statistical expectation of the gradient energy density within the physical integral scale.

[0087] The specific operation of eigenvalue decomposition to obtain structural tensor eigenvalues ​​is as follows: the structural tensor field is decomposed into three non-negative structural tensor eigenvalues ​​and three mutually orthogonal structural tensor eigenvectors. The structural tensor eigenvalue with the largest value among the three structural tensor eigenvalues ​​represents the structural strength along the normal direction of the anatomical structure, and the two structural tensor eigenvalues ​​with the smaller values ​​among the three structural tensor eigenvalues ​​represent the structural strength along the tangential direction of the anatomical structure. By comparing the numerical differences among the three structural tensor eigenvalues, flat regions, planar structural regions, and tubular structural regions in three-dimensional medical volumetric imaging data are distinguished.

[0088] This step uses the Jacobi iteration method to perform eigenvalue decomposition on the structure tensor field, obtaining three structure tensor eigenvalues ​​arranged in descending numerical order.

[0089] When the values ​​of the three structural tensor eigenvalues ​​are all below 10 square Henle units per square millimeter, it corresponds to a homogeneous region such as liver parenchyma; when the largest structural tensor eigenvalue is greater than 400 square Henle units per square millimeter and significantly greater than the other two structural tensor eigenvalues, it corresponds to a planar structural region such as organ capsule. This quantitative analysis based on the differences in eigenvalue values ​​provides clear geometric constraints for subsequent anisotropic diffusion.

[0090] Furthermore, in step S1, the specific operation of constructing the anisotropic diffusion impedance field based on the structural tensor eigenvalues ​​is as follows:

[0091] An anisotropic diffusion tensor is constructed, which is the anisotropic diffusion impedance field. The eigenvectors of the anisotropic diffusion tensor are set to be consistent with the eigenvectors of the structure tensor of the structure tensor field. The eigenvalues ​​of the anisotropic diffusion tensor are defined as heat flux conductivity.

[0092] By constructing a new 3x3 symmetric matrix as the anisotropic diffusion impedance field, the characteristic basis of the anisotropic diffusion impedance field is forced to completely coincide with the characteristic basis of the structural tensor field, ensuring that the diffusion process is strictly decomposed along the normal and tangential coordinate systems of the anatomical structure, rather than along the Cartesian coordinate system of the image acquisition.

[0093] A tissue contrast sensitivity threshold is set, which is used to distinguish random thermal noise from solid anatomical boundaries at the energy level.

[0094] In this embodiment, the specific value of the tissue contrast sensitivity threshold is set to 900 square Henle units per square millimeter (i.e., the square of 30 Henle units per millimeter). The basis for choosing 900 square Henle units per square millimeter as the tissue contrast sensitivity threshold is that in abdominal enhanced CT images, the gray-scale difference between weak boundaries between soft tissues is usually between 20 and 40 Henle units. Setting this threshold can ensure that the algorithm identifies fluctuations below this energy level as noise and abrupt changes above this energy level as anatomical boundaries.

[0095] Calculate the square of the difference between the largest and smallest eigenvalues ​​of the three structure tensors, and define the square of the difference as the coherence index.

[0096] By performing subtraction and exponentiation operations, the square of the difference between the eigenvalues ​​of the maximum and minimum structure tensors is calculated to obtain the coherence index. The coherence index mathematically amplifies the directional differences, making the signal intensity of the anatomical boundary with a clear direction significantly higher than the intensity of the isotropic noise signal.

[0097] For the heat flux conductivity along the anatomical structure normal, a negative exponential decay function based on the coherence index and the tissue contrast sensitivity threshold is used for calculation. This makes the heat flux conductivity along the anatomical structure normal approach zero when the coherence index is greater than the tissue contrast sensitivity threshold to form an adiabatic barrier, and makes the heat flux conductivity along the anatomical structure normal approach one when the coherence index is less than the tissue contrast sensitivity threshold to allow heat flow.

[0098] The heat flux conductivity along the anatomical structure normal is calculated using a negative exponential function model. The specific calculation logic is as follows: calculate the ratio of the coherence index to the sum of the tissue contrast sensitivity threshold and the numerical stability constant, take the negative number of this ratio as the exponent of the natural exponential function, and the calculated function value is the heat flux conductivity along the anatomical structure normal.

[0099] In this embodiment, the numerical stability constant is set to one part per million to prevent the denominator from being zero. When the coherence index reaches twice the tissue contrast sensitivity threshold, the heat flux conductivity along the anatomical structure normal drops to below 0.135, physically creating a strong blockage to normal diffusion, thereby protecting edge sharpness; when the coherence index is close to zero, the heat flux conductivity along the anatomical structure normal is close to 1, allowing for smooth noise reduction.

[0100] For the heat flux conductivity along the anatomical structure tangentially, the value of the heat flux conductivity along the anatomical structure tangentially is kept constant at 1 to allow omnidirectional diffusion.

[0101] The method sets the eigenvalues ​​corresponding to the two tangential eigenvectors in the anisotropic diffusion impedance field to a constant value of 1.

[0102] The specific operation of performing directional smoothing on 3D medical volumetric image data by substituting the anisotropic diffusion impedance field into the diffusion equation is as follows: establish a diffusion evolution equation that follows the principle of non-equilibrium thermodynamics, use the anisotropic diffusion impedance field to correct the diffusion flux in the diffusion evolution equation, and as physical time progresses, use natural diffusion to eliminate noise in the flat areas of the 3D medical volumetric image data, and use normal blocking to maintain the boundary in the edge areas of the 3D medical volumetric image data, until the diffusion evolution equation reaches a quasi-steady state, and output a preprocessed image volume containing an enhanced gradient vector field.

[0103] The diffusion evolution equation is solved iteratively using the finite difference method with a time step of 0.2. In each iteration, the product of the current anisotropic diffusion impedance field and the current gradient field is calculated to obtain the diffusion flux. Then, the divergence of the diffusion flux is calculated to update the voxel values ​​of the three-dimensional medical volumetric image data.

[0104] After twenty to fifty iterations, the diffusion evolution equation reaches a quasi-steady state. At this point, the output preprocessed image volume is more than 30% better than the original three-dimensional medical volume image data in terms of signal-to-noise ratio improvement by decibels compared with traditional Gaussian filtering, and the edge gradient preservation exponential is improved by 40%.

[0105] The output enhanced gradient vector field is strictly directed towards the normal of the anatomical structure. Furthermore, due to tangential diffusion at weak boundaries, the vector direction is made continuous and enhanced, providing a high signal-to-noise ratio data foundation for extracting an accurate initial seed point set in step S2 and constructing a high-precision gradient driving term in step S3.

[0106] Step S2: Receive the preprocessed image volume output from step S1, construct a Hessian matrix for the voxels in the preprocessed image volume, calculate the eigenvalues ​​of the Hessian matrix, remove tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix, and extract regions exhibiting grayscale energy convergence characteristics as the initial seed point set.

[0107] In this embodiment, after receiving the high signal-to-noise ratio preprocessed image volume output in step S1, step S2 is used to solve the problem of locating abnormal structures in an anatomical environment of unknown size. This step adopts differential geometry and linear scale space theory, and through the topological invariance of the analytical second derivative, accurately locks the three-dimensional convergence singularity of gradient energy in multi-scale space, providing an absolutely safe initial seed for subsequent segmentation.

[0108] Specifically, in step S2, the Hessian matrix is ​​constructed for the voxels in the preprocessed image volume, and the eigenvalues ​​of the Hessian matrix are calculated as follows:

[0109] A discrete physical observation aperture sequence containing multiple Gaussian smoothing scale parameters with different values ​​is defined, and a corresponding three-dimensional Gaussian kernel function is constructed for each Gaussian smoothing scale parameter in the discrete physical observation aperture sequence.

[0110] Given that the physical size of lesions (such as tumors and nodules) is an unknown random variable, in order to capture all possible lesions, the method must construct a series of observers with different physical apertures.

[0111] In this embodiment, the discrete physical observation aperture sequence is set as a set of values ​​starting from 1.0 mm, increasing in steps of 0.5 mm to 5.0 mm. This range is based on clinical statistical data and covers common diameter ranges from small lung nodules (usually less than 3 mm in diameter) to large liver tumors (usually more than 3 cm in diameter). For each value in the sequence, the method constructs a corresponding three-dimensional Gaussian kernel function, and the standard deviation of the three-dimensional Gaussian kernel function is the currently selected Gaussian smoothing scale parameter.

[0112] A three-dimensional convolution operation is performed on the volume of the preprocessed image using a three-dimensional Gaussian kernel function to generate a multi-scale spatial response field. The square of the Gaussian smoothing scale parameter is calculated and defined as the scale normalization factor. The multi-scale spatial response field is then subjected to a multiplicative weighted operation using the scale normalization factor to generate a scale-normalized multi-scale spatial response field. Finally, the second-order partial derivative of each voxel in the scale-normalized multi-scale spatial response field in the Cartesian coordinate system is calculated.

[0113] Specifically, this step involves convolving a three-dimensional Gaussian kernel function with the volume of the preprocessed image. This process is physically equivalent to smoothing the image using apertures of different sizes.

[0114] It is worth noting that as the Gaussian smoothing scale parameter increases, the magnitude of the second derivative of the image will decrease quadratically, resulting in a weaker response signal for large-scale lesions. To address this physical distortion, the method introduces a scale normalization factor, which is the square of the Gaussian smoothing scale parameter (corresponding to a scale normalization exponent of two).

[0115] The method multiplies the multi-scale spatial response field by a scale normalization factor, which serves as an energy compensation in the algorithm.

[0116] A 3×3 second-order symmetric matrix describing the curvature of the local energy surface is constructed using second-order partial derivatives. The 3×3 second-order symmetric matrix is ​​defined as the Hessian matrix. Eigenvalue decomposition is performed on the Hessian matrix to obtain three Hessian matrix eigenvalues ​​representing the principal curvature. The three Hessian matrix eigenvalues ​​are then arranged in ascending order according to their absolute values.

[0117] This step calculates the second-order partial derivatives and mixed partial derivatives of each voxel with respect to the three axes of the Cartesian coordinate system, and combines them into a 3×3 symmetric matrix, namely the Hessian matrix.

[0118] The Hessian matrix physically describes the curvature of a local energy surface in three-dimensional space, i.e., the second-order Taylor expansion coefficients. Subsequently, the method uses a linear algebra algorithm to perform eigenvalue decomposition on the Hessian matrix, obtaining three eigenvalues. In order to unify the physical meaning and eliminate ambiguity, the method sorts these three eigenvalues ​​in ascending order of their absolute values.

[0119] The Hessian matrix eigenvalue with the smallest absolute value among the three Hessian matrix eigenvalues ​​is defined as the first principal curvature along the structural axis, the Hessian matrix eigenvalue with the largest absolute value among the three Hessian matrix eigenvalues ​​is defined as the third principal curvature along the cross-sectional direction, and the Hessian matrix eigenvalue with the middle absolute value among the three Hessian matrix eigenvalues ​​is defined as the second principal curvature.

[0120] This sequencing step is crucial because the first principal curvature corresponds to the direction of slowest change in local geometry (i.e., along the axis of the tubular structure), while the third principal curvature corresponds to the direction of fastest change (i.e., the cross-sectional direction). This explicit physical definition provides a mathematical basis for the subsequent removal of blood vessels and bone walls.

[0121] Furthermore, in step S2, the specific operation of removing tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix and extracting regions exhibiting gray-level energy convergence characteristics as the initial seed point set is as follows:

[0122] Calculate the absolute value of the first principal curvature and define the absolute value of the first principal curvature as the minimum curvature intensity term. The minimum curvature intensity term is used to ensure that the voxel has significant curvature in all directions to eliminate planar structures.

[0123] In physical geometry, for planar structures (such as organ capsules), the first principal curvature approaches zero; for tubular structures (such as blood vessels), the first principal curvature also approaches zero (no curvature along the axial direction of the blood vessel). Only when the absolute value of the first principal curvature is large enough does it mean that the voxel has undergone significant curvature in all three directions, that is, the voxel is located at the core of a spherical structure. Therefore, this step uses the absolute value of the first principal curvature as the criterion for distinguishing between a true sphere and a blood vessel section or planar artifact.

[0124] The ratio of the absolute value of the first principal curvature to the absolute value of the third principal curvature is calculated. This ratio is defined as the geometric isotropic factor, which is used to quantify the degree to which the local geometric structure approximates a sphere in order to eliminate tubular structures.

[0125] The physical meaning of the geometric isotropic factor is to describe the uniformity of local curvature. For an ideal sphere, the three principal curvatures are equal, and the geometric isotropic factor is 1. For a tubular structure, the first principal curvature is much smaller than the third principal curvature, and the geometric isotropic factor approaches zero.

[0126] In this embodiment, in order to prevent division by zero errors caused by the third principal curvature approaching zero in flat regions, the method introduces a very small numerical stability constant in the denominator, with a value of one part per million.

[0127] Meanwhile, to adjust the tolerance for irregular lesion shapes, this invention also introduces a morphological selectivity index as a power, which is preferably set between 0.5 and 2.0 in this embodiment. The larger the morphological selectivity index, the more stringent the algorithm's requirement for perfect circles.

[0128] Define a polarity-gated operator, which is used to detect the sign consistency of the three Hessian matrix eigenvalues. When it is necessary to extract bright lesions, configure the polarity-gated operator to output a non-zero value when all three Hessian matrix eigenvalues ​​are negative. When it is necessary to extract dark lesions, configure the polarity-gated operator to output a non-zero value when all three Hessian matrix eigenvalues ​​are positive.

[0129] The polarity gating operator endows the algorithm with the ability to adapt to the polarity of lesions. For bright tumors in enhanced CT, the local gray level shows a peak shape, and all three principal curvatures are negative. For low-density cysts in plain CT, the gray level shows a trough shape, and all three principal curvatures are positive. By configuring the polarity gating operator, it is possible to flexibly switch between bright target capture and dark target capture modes, which solves the problem of detecting heterogeneous lesions. This is something that cannot be achieved by existing single threshold segmentation methods.

[0130] The product of the minimum curvature intensity term, the geometric isotropic factor, and the polarity gate operator is calculated, and the product is defined as the structural significance response value of the voxel under the Gaussian smoothing scale parameter.

[0131] The method performs a multiplication operation, merging the three physical factors into a single scalar value. The structural saliency response value comprehensively reflects the probability that the voxel belongs to the core of the lesion at the current scale.

[0132] The maximum value of the structural significance response of the voxel under all Gaussian smoothing scale parameters is selected and defined as the multi-scale convergence probability field, which represents the probability that the voxel belongs to the core of the lesion.

[0133] By traversing all scales in the discrete physical observation aperture sequence and taking the maximum response value of each voxel at different scales, this maximum likelihood response strategy ensures that the strongest signal response can be obtained at the best matching scale regardless of the physical size of the lesion, thereby achieving robust detection of lesions at multiple scales.

[0134] Calculate the non-zero mean and non-zero standard deviation of the multi-scale convergence probability field, and calculate the sum of the products of the non-zero mean, non-zero standard deviation, and confidence factor. Define the sum as the adaptive truncation threshold. Extract the regions in the multi-scale convergence probability field whose values ​​are greater than the adaptive truncation threshold, generate a binarized mask, and define the binarized mask as the initial seed point set.

[0135] To eliminate low-probability background noise, this invention performs a statistical truncation operation. In this embodiment, the confidence factor is preferably set between 2.0 and 3.0. This value is based on the statistical regularity derived from the normal distribution and aims to extract high-confidence regions with a significance level of over 95%. The method marks voxels in the multi-scale convergent probability field that are greater than the adaptive truncation threshold as one, and the rest as zero, generating a binarized mask.

[0136] The spatial region corresponding to the binarized mask is the initial seed point set. The initial seed point set avoids edges and adhesion areas and is accurately located at the topological core of the lesion. Compared with the traditional manual interactive seed point selection method, the method of generating the initial seed point set in this embodiment not only shortens the processing time from an average of five minutes to less than three seconds, but also eliminates the operator's subjective bias, providing stable and high-confidence initialization conditions for the automated segmentation based on the level set evolution equation in step S3.

[0137] Step S3: Receive the initial seed point set output in step S2, construct an initial level set function with the initial seed point set as the zero level set region, construct a curvature inverse braking term, construct a gradient driving term based on the enhanced gradient vector field in step S1, establish a level set evolution equation composed of the curvature inverse braking term and the gradient driving term, perform iterative evolution on the initial level set function using the level set evolution equation, stop local evolution when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function makes the curvature inverse braking term greater than the gradient driving term, and output the evolved level set function.

[0138] In this embodiment, after receiving the high-confidence initial seed point set output by step S2, step S3 uses the enhanced gradient vector field output by step S1 as an external potential energy guide, aiming to solve the technical problems of leakage of three-dimensional anatomical structures in weak boundary regions and difficulty in accurately capturing topological structures.

[0139] Step S3, based on fluid mechanics and variational geometric flow theory, simulates the dynamic deformation process of a controlled pressurized elastic film in a non-uniform medium. The evolution process is manifested as a dynamic game of three physical forces: vector advection traction force, scalar expansion thrust, and Yang-Laplace curvature reverse braking force.

[0140] Specifically, in step S3, the initial level set function is constructed with the initial seed point set as the zero level set region, and the curvature inverse braking term is constructed. The specific operations for constructing the gradient driving term based on the enhanced gradient vector field in step S1 are as follows:

[0141] Calculate the shortest Euclidean distance from each voxel in 3D space to the edge of the initial seed point set. Construct a symbolic distance field based on the shortest Euclidean distance. The symbolic distance field is defined to take negative values ​​inside the initial seed point set, positive values ​​outside the initial seed point set, and zero values ​​at the edge of the initial seed point set. Define the symbolic distance field as the initial level set function.

[0142] As a high-dimensional implicit scalar field, the initial level set function can handle the splitting and merging of topological structures. In this embodiment, the fast-progression algorithm is used to calculate the shortest Euclidean distance to ensure that the gradient magnitude of the initial level set function is strictly equal to 1. The symbolic distance field is used instead of the binary mask as the evolution object, which solves the numerical dissipation and jaggedness problems caused by the binary mask during the evolution process, so that the evolution interface maintains sub-voxel level smoothness.

[0143] Define a narrow band region where the initial level set function is close to zero. Update the initial level set function only within the narrow band region. Calculate the vector magnitude of the enhanced gradient vector field output in step S1. Construct an edge indicator scalar function based on the vector magnitude of the enhanced gradient vector field. Set the edge indicator scalar function to approach zero at strong edges where the vector magnitude of the enhanced gradient vector field is large, and set the edge indicator scalar function to approach 1 at flat regions where the vector magnitude of the enhanced gradient vector field is small.

[0144] In this embodiment, the width of the narrowband region is set to a range of three to five voxels. The basis for selecting the width of the narrowband region is that if the width of the narrowband region is less than three voxels, the evolution interface may cross the boundary of the computational domain in a single iteration, resulting in numerical truncation error; if the width of the narrowband region is greater than five voxels, the computational load will increase cubically.

[0145] For the edge indicator scalar function, the method uses an inverse square function for construction, which introduces a gradient cutoff threshold parameter. In this embodiment, the gradient cutoff threshold is set to 30 Henle units per millimeter. The selection of this value is based on the gradient statistical distribution of soft tissue and background noise: when the vector magnitude of the locally enhanced gradient vector field is significantly greater than 30 Henle units per millimeter, the value of the edge indicator scalar function drops sharply to below 0.1, which physically manifests as an impermeable barrier layer; otherwise, it manifests as a freely diffusing medium.

[0146] The specific operation of constructing the gradient driving term based on the enhanced gradient vector field in step S1 is as follows: the gradient driving term is decoupled into two independent parts: a vector advection traction term and a scalar expansion thrust term. The inner product of the enhanced gradient vector field and the gradient vector of the initial level set function is calculated, and the inner product is defined as the vector advection traction term. The vector advection traction term drives the zero level set interface to adhere to the anatomical structure boundary along the direction of the enhanced gradient vector field. The product of the preset expansion pressure coefficient and the edge indicator scalar function is calculated, and the product is defined as the scalar expansion thrust term. The scalar expansion thrust term simulates the fluid injection pressure driving the zero level set interface to expand outward.

[0147] The vector advection traction term utilizes the directional information of the enhanced gradient vector field to exert a traction effect on the zero-level set interface, forcing the normal direction of the zero-level set interface to be consistent with the edge normal direction of the anatomical structure, thereby solving the edge slip problem in traditional scalar evolution. The scalar expansion thrust term provides continuous outward power.

[0148] In this embodiment, the expansion pressure coefficient is set to a dimensionless constant between 1.0 and 5.0. The setting of the expansion pressure coefficient is based on the balance between Courant-Friedrich-Levy stability conditions and evolutionary efficiency: if the expansion pressure coefficient is too small, the evolution will stagnate at a local minimum; if the expansion pressure coefficient is too large, it will destroy numerical stability. By setting this range, it is ensured that the zero-level set interface can advance outward at a constant speed within the flat anatomical structure.

[0149] The specific operation of constructing the curvature inverse braking term is as follows: Based on Yang Laplace's law, an adaptive surface tension mechanism is constructed. The normalized divergence of the gradient vector of the initial level set function is calculated. The normalized divergence of the gradient vector is defined as the local mean curvature of the zero level set interface. A surface tension coefficient is introduced and set as a function of the edge indicator scalar function. This allows the surface tension coefficient to automatically increase at weak boundaries where the edge indicator scalar function has a large value. The product of the surface tension coefficient, the edge indicator scalar function, and the local mean curvature is calculated and defined as the curvature inverse braking term.

[0150] According to Laplace's law, the additional pressure on a curved liquid surface is proportional to the local average curvature. This invention utilizes this physical law to construct an adaptive mechanism: in flat regions, the surface tension coefficient is maintained at a base value of 0.1, allowing the interface to deform freely; while in adhesion regions with weak gradients (where the edge indicator scalar function value is close to one), the surface tension coefficient automatically increases to five to ten times the base value. Comparative data shows that after introducing the adaptive surface tension mechanism, the leakage rate in weak boundary regions such as blood vessel bifurcation is reduced by 92%. This dynamic adjustment mechanism enhances the anti-leakage capability as the gradient weakens, physically forming a dynamic blockage of narrow leakage channels.

[0151] Furthermore, in step S3, an evolution equation for the level set, composed of a curvature inverse braking term and a gradient driving term, is established. The initial level set function is iteratively evolved using this equation. Local evolution stops when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function causes the curvature inverse braking term to be greater than the gradient driving term. The specific operation for outputting the evolved level set function is as follows:

[0152] A partial differential evolution equation following the Hamiltonian-Jacobi form is established. The partial differential evolution equation describes the dynamic deformation process of the initial level set function with physical time. A regularized Dirac function is introduced into the partial differential evolution equation to restrict the evolutionary force to the physical location of the zero level set interface.

[0153] The regularized Dirac function strictly restricts all physical forces to a voxel range near the zero-level set interface, ensuring the locality of energy calculation and avoiding interference from background noise far from the interface in the evolution process.

[0154] The evolution rate of the partial differential evolution equation is determined by the algebraic sum of the vector advection traction term, the scalar expansion thrust term, and the curvature reverse braking term. During the iterative evolution process, for the flat internal region of the anatomical structure, the scalar expansion thrust term dominates, driving the zero level set interface to grow outward.

[0155] In the flat interior region, due to the weak gradient, the values ​​of the vector advection traction term and the curvature reverse braking term are close to zero. The method manifests as simple volume expansion, achieving rapid filling of the lesion interior.

[0156] For the strong edge regions of the anatomical structure, the edge indicator scalar function value approaches zero, causing the evolution rate to approach zero, locking the zero-level set interface on the boundary of the anatomical structure. For the weak boundary adhesion regions of the anatomical structure, the narrow channel causes the local average curvature value to increase sharply, so that the contraction force generated by the curvature reverse braking term exceeds the expansion force generated by the scalar expansion thrust term, causing the local evolution rate to become negative, forcing the zero-level set interface to automatically shrink back and stop growing.

[0157] This is the leakage prevention mechanism of this embodiment. When the zero-level set interface attempts to pass through a narrow vascular adhesion or tissue gap, the local geometry of the zero-level set interface will bend violently, causing the local average curvature value to increase sharply. The value of the curvature reverse braking term will increase exponentially, instantly exceeding the constant scalar expansion thrust term. At this time, the total evolution rate calculated by the partial differential evolution equation changes from positive to negative. Physically, this is manifested as the interface being subjected to huge inward surface tension and undergoing elastic retraction, thereby achieving automatic truncation at the weak boundary.

[0158] The spatial derivative of the initial level set function is discretized using the central difference method, and the temporal derivative of the initial level set function is discretized using the forward Euler method. The time step is dynamically adjusted according to the numerical stability condition, which requires that the time step be less than the minimum value of the ratio of the spatial grid size to the maximum evolution rate.

[0159] The numerical stability condition is the Courant-Friedrich-Levy condition. In this embodiment, the time step is not a fixed value, but is dynamically determined by the maximum flow velocity of the entire field in the current iteration, and is usually limited to below 0.25. This constraint ensures that the numerical information will not cross more than one grid cell in a single iteration, thereby guaranteeing the convergence and stability of the partial differential evolution equation solution and avoiding numerical oscillations and interface breakage caused by excessive time steps.

[0160] The iterative evolution is continuously executed until the partial differential evolution equation reaches a steady state or the total energy change rate is lower than a preset threshold. The final steady-state scalar field is output as the evolved level set function, and the zero level set interface of the evolved level set function is extracted as the surface of the segmented three-dimensional anatomical structure.

[0161] By monitoring the total energy change rate of the initial level set function in five consecutive iterations, when the change rate is less than one-thousandth, the system is determined to have reached a steady state of mechanical equilibrium. At this point, the evolved level set function output contains complete geometric and topological information of the three-dimensional anatomical structure.

[0162] By extracting the zero level set interface of the evolved level set function, a high-precision three-dimensional surface is obtained after segmentation. While ensuring sub-voxel level accuracy, this three-dimensional surface perfectly eliminates tubular adhesions and sheet artifacts by using the curvature inverse braking term, realizing topological adaptive segmentation of complex anatomical structures, and providing topologically correct and geometrically accurate input data for the mesh reconstruction based on the physical model in step S4.

[0163] Step S4: Receive the evolved level set function output from step S3, extract the initial mesh vertices and initial vertex connection relationships from the zero level set interface of the evolved level set function, calculate the data attachment energy term for the initial mesh vertices based on the enhanced gradient vector field in step S1, calculate the elastic deformation energy term for the initial mesh vertices based on the initial vertex connection relationships, construct a variational energy functional composed of the weighted sum of the data attachment energy term and the elastic deformation energy term, calculate the optimized spatial coordinates of the initial mesh vertices by solving the minimum value of the variational energy functional, and generate a three-dimensional mesh model based on the optimized spatial coordinates.

[0164] In this embodiment, after receiving the level set function after steady-state evolution output from step S3, step S4 aims to transform the implicit topology into a high-precision explicit 3D mesh model. Step S4 extracts discrete explicit geometric representations from the continuous implicit field. To address the staircase artifacts and volume collapse caused by mesh smoothing in traditional isosurface extraction algorithms, this embodiment constructs a physical model based on adaptive generalized potential functionals. This step treats the discrete isosurface mesh as a point mass system placed in a non-uniform elastic medium, utilizing the dynamic balance of three physical forces—gradient adsorption force, tangential tension, and adaptive normal pressure—to achieve sub-voxel-level reconstruction of the 3D anatomical structure.

[0165] Specifically, in step S4, the evolved level set function output from step S3 is received, and the initial mesh vertices and initial vertex connection relationships are extracted from the zero level set interface of the evolved level set function as follows:

[0166] The zero level set interface in the evolved level set function output by step S3 is parsed using the moving cube isosurface extraction algorithm. The zero level set interface is parsed to obtain the discrete initial mesh vertex set and the initial vertex connection relationship. The initial mesh vertex set contains the initial mesh vertex coordinate information defined in three-dimensional space, and the initial vertex connection relationship contains the triangle mesh list that defines the physical topological adjacency attributes between the initial mesh vertices.

[0167] In this embodiment, the method uses the moving cube algorithm as a topological discretization tool. The method first divides the three-dimensional space into a continuous cubic voxel grid. The resolution of the voxel grid is consistent with the physical resolution of the original three-dimensional medical volumetric image data. Then, the method detects the signed distance value of the eight corner points of each voxel. When the signed distance value of the voxel corner point crosses zero, the method determines that the zero level set interface passes through the voxel and calculates the specific coordinate position of the initial grid vertex on the voxel edge according to linear interpolation. At this time, although the initial grid vertex set generated has the correct topological structure, due to the limitations of the discrete sampling characteristics of the original image (e.g., the layer thickness is 0.625 mm), the surface formed by the initial grid vertex set often exhibits non-physical step-like artifacts and is in a high potential physical excited state. A subsequent variational energy relaxation process is needed to release stress and improve geometric accuracy.

[0168] The specific operation for calculating the data-attached energy term at the initial grid vertices based on step S1 is as follows: For the enhanced gradient vector field defined on the discrete voxel grid points by the output of step S1, a trilinear interpolation operator is constructed. The trilinear interpolation operator is used to establish a mapping relationship from the discrete data space to the continuous physical space, and the discrete enhanced gradient vector field is transformed into a continuous gradient potential energy density field.

[0169] The trilinear interpolation operator constructed by this method can calculate the gradient vector at any continuous coordinate position in 3D space by performing a distance-weighted average of the gradient vectors of the eight adjacent voxel grid points surrounding the target point. This operation mathematically constructs a continuous vector field with infinite resolution, providing a data foundation for the initial grid vertices to break through the discrete limitations of voxel grids. Comparative experiments show that if only nearest neighbor interpolation is used, the reconstruction error is locked at 0.5 voxels; while after using the trilinear interpolation operator, the average geometric error of the reconstructed surface is reduced to 0.08 voxels, proving the necessity of continuous field construction for sub-voxel level reconstruction.

[0170] For each initial mesh vertex, the interpolation gradient vector of the initial mesh vertex at continuous physical space coordinates is calculated using the trilinear interpolation operator. The square of the magnitude of the interpolation gradient vector is calculated, and the negative value of the square of the magnitude of the interpolation gradient vector is defined as the data attachment energy term. The data attachment energy term physically constructs a continuous gravitational potential well. The energy valley of the gravitational potential well strictly corresponds to the gradient maximum position of the real physical boundary of the anatomical structure. The gravitational potential well generates a gradient attraction force pointing to the gradient maximum position. The gradient attraction force drives the initial mesh vertex to break through the voxel resolution limitation and automatically slides the initial mesh vertex to the position of local gradient energy maximum in the continuous physical space with sub-voxel precision.

[0171] The data attachment energy term essentially defines an external potential energy field in which each initial grid vertex is treated as a force-bearing point mass, and the gradient adsorption force it experiences is proportional to the negative gradient of the data attachment energy term. Physically, this means that the initial grid vertices will automatically migrate to the local maximum location of the enhanced gradient vector field (i.e., the true edge of the anatomical structure). Due to the introduction of the trilinear interpolation operator, this migration is no longer limited to discrete voxel grid points, but can be fine-tuned in a continuous space at the sub-voxel level, thereby achieving high-precision reconstruction that surpasses the resolution of the original CT image.

[0172] Further, in step S4, the elastic deformation energy term for the initial mesh vertices is calculated based on the initial vertex connection relationship. A variational energy functional is constructed, consisting of a weighted sum of the data attachment energy term and the elastic deformation energy term. The optimized spatial coordinates of the initial mesh vertices are calculated by solving for the minimum value of the variational energy functional. The specific operation of generating a 3D mesh model based on the optimized spatial coordinates is as follows:

[0173] The mean curvature vector of the initial mesh vertices is calculated using the discrete Laplace Beltrami operator. An area-weighted normal vector is introduced to prevent numerical singularities when the initial mesh vertices degenerate. The mean curvature vector is orthogonally decomposed into a tangential relaxation component and a normal smoothing component.

[0174] This invention calculates the discrete Laplace Beltrami operator using the cotangent weight formula. This operator can accurately describe the local geometric curvature of the triangular mesh surface. To prevent numerical calculation collapse when the area of ​​the triangular mesh approaches zero, an area weighting mechanism is introduced to stabilize the calculation of the normal vector. Subsequently, based on the orthogonal decomposition principle of forces, the average curvature vector is decomposed into a tangential relaxation component along the tangential plane direction and a normal smoothing component along the normal direction, providing a physical basis for subsequent independent control.

[0175] Tangential tension is constructed based on tangential relaxation components. Tangential tension optimizes the shape quality and vertex distribution uniformity of the triangular mesh in the initial vertex connection relationship. Tangential tension does not change the volume of the 3D mesh model.

[0176] In layman's terms, the effect of tangential tension is similar to covering the mesh surface with a frictionless elastic membrane. Tangential tension drives the initial mesh vertices to slide on the tangential plane, causing the triangular mesh to tend towards an equilateral triangle, thereby eliminating low-mass elements such as narrow triangles. Since tangential motion does not produce displacement components perpendicular to the surface, this process strictly follows the law of volume conservation, avoiding the model volume shrinkage problem caused by traditional smoothing algorithms.

[0177] Normal pressure is constructed based on normal smoothing components, which eliminates high-frequency noise on the initial mesh vertex surface. Local adaptive stiffness coefficients are constructed, and the local adaptive stiffness coefficients are set to be inversely proportional to the local magnitude of the enhanced gradient vector field. The local adaptive stiffness coefficients are set to approach zero at strong edges where the local magnitude of the enhanced gradient vector field is large in order to turn off the normal pressure, and the local adaptive stiffness coefficients are set to increase at weak edges where the local magnitude of the enhanced gradient vector field is small in order to maintain the shape using surface tension.

[0178] This is the mechanism by which this embodiment resolves the contradiction between overfitting at strong edges and volume collapse at weak edges. The local adaptive stiffness coefficient constructed in this invention simulates the physical properties of a non-uniform medium.

[0179] In this embodiment, the local adaptive stiffness coefficient is set to a range of 0.1 to 1.0, based on the following criteria: When the magnitude of the local enhanced gradient vector field is greater than 50 Henle units per millimeter (strong edge), the local adaptive stiffness coefficient is forcibly set to below 0.1, causing the normal pressure to fail and the initial mesh vertices to be completely controlled by the gradient adsorption force, thereby preserving sharp anatomical features; when the magnitude of the local enhanced gradient vector field is less than 10 Henle units per millimeter (weak edge or flat region), the local adaptive stiffness coefficient is linearly increased to 1.0, causing the normal pressure to dominate, using surface tension to eliminate step artifacts and maintain a smooth geometry.

[0180] The algebraic sum of the tangential tension and the normal pressure weighted by the local adaptive stiffness coefficient is calculated and defined as the elastic deformation energy term. A variational energy functional containing the data attachment energy term and the elastic deformation energy term is constructed. The variational energy functional is iteratively solved using the semi-implicit Euler method until the maximum displacement of the initial mesh vertex is less than the preset thermodynamic convergence threshold.

[0181] In this embodiment, the thermodynamic convergence threshold is set to 0.01 mm. In each iteration, the displacement of all initial mesh vertices is calculated. When the maximum displacement is less than the threshold, the system is determined to have reached mechanical equilibrium (i.e., frozen state). The use of semi-implicit Euler method ensures numerical stability at large time steps (e.g., 0.2), which significantly improves the solution efficiency. It can usually converge within twenty iterations.

[0182] The output is the optimized spatial coordinates of the initial mesh vertex when it reaches mechanical equilibrium. The initial vertex connection relationship remains unchanged. The coordinates of the initial mesh vertex are updated to the optimized spatial coordinates to generate a 3D mesh model with sub-voxel level geometric accuracy and smooth topology.

[0183] The final generated 3D mesh model not only has a geometric error of less than 0.1 mm (compared to the gold standard), but also has a high-quality triangular mesh topology (no self-intersection, no degradation). Through this step, the present invention solves the problem of crossing from discrete voxel data to a continuous physical model, eliminates the limitation of image resolution on model accuracy, and the output mesh model can be directly used for subsequent finite element biomechanical analysis or high-precision 3D printing, providing a reliable physical model basis for clinical surgical planning.

[0184] Example 2

[0185] like Figure 2 As shown, the present invention also discloses a segmentation and reconstruction system for three-dimensional medical image data, comprising: an image preprocessing module, a seed point extraction module, a level set evolution module, and a mesh model generation module, wherein;

[0186] The image preprocessing module receives three-dimensional medical volumetric image data, calculates the gradient information of voxels in the three-dimensional medical volumetric image data to construct a structure tensor field, performs eigenvalue decomposition on the structure tensor field to obtain structure tensor eigenvalues, constructs an anisotropic diffusion impedance field based on the structure tensor eigenvalues, substitutes the anisotropic diffusion impedance field into the diffusion equation to perform directional smoothing processing on the three-dimensional medical volumetric image data, and outputs a preprocessed image volume containing an enhanced gradient vector field.

[0187] The seed point extraction module receives the preprocessed image volume output by the image preprocessing module, constructs a Hessian matrix for the voxels in the preprocessed image volume, calculates the eigenvalues ​​of the Hessian matrix, removes tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix, and extracts regions exhibiting gray-scale energy convergence characteristics as the initial seed point set.

[0188] The level set evolution module receives the initial seed point set output by the seed point extraction module, constructs an initial level set function with the initial seed point set as the zero level set region, constructs a curvature inverse braking term, constructs a gradient driving term based on the gradient vector field enhanced by the image preprocessing module, establishes a level set evolution equation composed of the curvature inverse braking term and the gradient driving term, performs iterative evolution on the initial level set function using the level set evolution equation, stops local evolution when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function makes the curvature inverse braking term greater than the gradient driving term, and outputs the evolved level set function.

[0189] The mesh model generation module receives the evolved level set function output by the level set evolution module, extracts the initial mesh vertices and initial vertex connections from the zero level set interface of the evolved level set function, calculates the data attachment energy term for the initial mesh vertices based on the enhanced gradient vector field of the image preprocessing module, calculates the elastic deformation energy term for the initial mesh vertices based on the initial vertex connections, constructs a variational energy functional composed of the weighted sum of the data attachment energy term and the elastic deformation energy term, calculates the optimized spatial coordinates of the initial mesh vertices by solving the minimum value of the variational energy functional, and generates a three-dimensional mesh model based on the optimized spatial coordinates.

[0190] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for segmenting and reconstructing three-dimensional medical image data, characterized in that, Includes the following steps: Step S1: Receive three-dimensional medical volumetric image data, calculate the gradient information of voxels in the three-dimensional medical volumetric image data to construct a structure tensor field, perform eigenvalue decomposition on the structure tensor field to obtain structure tensor eigenvalues, construct an anisotropic diffusion impedance field based on the structure tensor eigenvalues, substitute the anisotropic diffusion impedance field into the diffusion equation to perform directional smoothing processing on the three-dimensional medical volumetric image data, and output a preprocessed image volume containing an enhanced gradient vector field. Step S2: Receive the preprocessed image volume output in step S1, construct a Hessian matrix for the voxels in the preprocessed image volume, calculate the eigenvalues ​​of the Hessian matrix, remove tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix, and extract regions exhibiting gray-scale energy convergence characteristics as the initial seed point set. Step S3: Receive the initial seed point set output in step S2, construct an initial level set function with the initial seed point set as the zero level set region, construct a curvature inverse braking term, construct a gradient driving term based on the enhanced gradient vector field in step S1, establish a level set evolution equation composed of the curvature inverse braking term and the gradient driving term, perform iterative evolution on the initial level set function using the level set evolution equation, stop local evolution when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function makes the curvature inverse braking term greater than the gradient driving term, and output the evolved level set function. Specifically, the gradient driving term is decoupled into two independent parts: a vector advection traction term and a scalar expansion thrust term. The inner product of the enhanced gradient vector field and the gradient vector of the initial level set function is calculated, and the inner product is defined as the vector advection traction term. The product of the preset expansion pressure coefficient and the edge indicator scalar function is calculated, and the product is defined as the scalar expansion thrust term. An adaptive surface tension mechanism is constructed based on Yang-Laplace's law. The normalized divergence of the gradient vector of the initial level set function is calculated. The normalized divergence of the gradient vector is defined as the local mean curvature of the zero level set interface. A surface tension coefficient is introduced and set as a function of the edge indicator scalar function. This allows the surface tension coefficient to automatically increase at weak boundaries where the edge indicator scalar function has a large value. The product of the surface tension coefficient, the edge indicator scalar function, and the local mean curvature is calculated and defined as the curvature inverse braking term. A partial differential evolution equation following the Hamiltonian Jacobian form is established, and the evolution rate of the partial differential evolution equation is determined by the algebraic sum of the vector advection traction term, the scalar expansion thrust term, and the curvature reverse braking term. Step S4: Receive the evolved level set function output from step S3. Extract the initial mesh vertices and initial vertex connections from the zero level set interface of the evolved level set function. Calculate the data attachment energy term for the initial mesh vertices based on the enhanced gradient vector field from step S1. Calculate the elastic deformation energy term for the initial mesh vertices based on the initial vertex connections. Construct a variational energy functional consisting of a weighted sum of the data attachment energy term and the elastic deformation energy term. Calculate the optimized spatial coordinates of the initial mesh vertices by solving for the minimum value of the variational energy functional. Generate a 3D mesh model based on the optimized spatial coordinates. Specifically, for each initial grid vertex, the interpolation gradient vector of the initial grid vertex at continuous physical space coordinates is calculated using the trilinear interpolation operator; the negative value of the square of the magnitude of the interpolation gradient vector is defined as the data attachment energy term; The mean curvature vector of the initial mesh vertices is calculated using the discrete Laplace Beltrami operator. The mean curvature vector is orthogonally decomposed into a tangential relaxation component and a normal smoothing component. Tangential tension is constructed based on the tangential relaxation component, and normal pressure is constructed based on the normal smoothing component. A local adaptive stiffness coefficient is constructed, which is set to be inversely proportional to the local modulus of the enhanced gradient vector field. The local adaptive stiffness coefficient is set to approach zero at strong edges with large local modulus of the enhanced gradient vector field to close the normal pressure, and is set to increase at weak edges with small local modulus of the enhanced gradient vector field to maintain the shape using surface tension. The algebraic sum of the tangential tension and the normal pressure weighted by the local adaptive stiffness coefficient is calculated, and the algebraic sum is defined as the elastic deformation energy term.

2. The method for segmenting and reconstructing three-dimensional medical image data according to claim 1, characterized in that, In step S1, the specific operation of calculating the gradient information of voxels in the three-dimensional medical volumetric image data to construct the structure tensor field is as follows: A physical observation scale is defined, which is used to filter out quantum noise in 3D medical volumetric image data. A Gaussian kernel function that conforms to the physical observation scale is used to perform 3D convolution smoothing on the 3D medical volumetric image data to generate a continuous differentiable physical field. Calculate the gradient vector of each voxel in a continuously differentiable physical field in three-dimensional space. The gradient vector indicates the direction of gray-scale change of the continuously differentiable physical field at the physical observation scale. Calculate the outer product matrix of the gradient vector and its transpose. The outer product matrix describes the tensor product of the gradient directions at a single point. Define a physical integral scale, which determines the range of the physical support domain of local geometric perception. Use a Gaussian weighting function that conforms to the physical integral scale to perform a weighted average calculation on the outer product matrix within the range of the physical support domain to generate a second-order symmetric positive semi-definite matrix describing the statistical distribution of local gradient energy. The second-order symmetric positive semi-definite matrix is ​​the structure tensor field. The specific operation of eigenvalue decomposition to obtain structural tensor eigenvalues ​​is as follows: the structural tensor field is decomposed into three non-negative structural tensor eigenvalues ​​and three mutually orthogonal structural tensor eigenvectors. The structural tensor eigenvalue with the largest value among the three structural tensor eigenvalues ​​represents the structural strength along the normal direction of the anatomical structure, and the two structural tensor eigenvalues ​​with the smaller values ​​among the three structural tensor eigenvalues ​​represent the structural strength along the tangential direction of the anatomical structure. By comparing the numerical differences among the three structural tensor eigenvalues, flat regions, planar structural regions, and tubular structural regions in three-dimensional medical volumetric imaging data are distinguished.

3. The method for segmenting and reconstructing three-dimensional medical image data according to claim 2, characterized in that, In step S1, the specific operation of constructing the anisotropic diffusion impedance field based on the structural tensor eigenvalues ​​is as follows: An anisotropic diffusion tensor is constructed, which is the anisotropic diffusion impedance field. The eigenvectors of the anisotropic diffusion tensor are set to be consistent with the eigenvectors of the structure tensor of the structure tensor field. The eigenvalues ​​of the anisotropic diffusion tensor are defined as heat flux conductivity. Set a tissue contrast sensitivity threshold, which is used to distinguish random thermal noise from solid anatomical boundaries at the energy level. Calculate the square of the difference between the largest and smallest structure tensor eigenvalues ​​among the three structure tensor eigenvalues, and define the square of the difference as the coherence index. For the heat flux conductivity along the anatomical structure normal, a negative exponential decay function based on the coherence index and the tissue contrast sensitivity threshold is used for calculation. This makes the heat flux conductivity along the anatomical structure normal approach zero when the coherence index is greater than the tissue contrast sensitivity threshold to form an adiabatic barrier, and makes the heat flux conductivity along the anatomical structure normal approach one when the coherence index is less than the tissue contrast sensitivity threshold to allow heat flow. For the heat flux conductivity along the tangential direction of the anatomical structure, the value of the heat flux conductivity along the tangential direction of the anatomical structure is kept constant at one to allow omnidirectional diffusion; The specific operation of performing directional smoothing on 3D medical volumetric image data by substituting the anisotropic diffusion impedance field into the diffusion equation is as follows: establish a diffusion evolution equation that follows the principle of non-equilibrium thermodynamics, use the anisotropic diffusion impedance field to correct the diffusion flux in the diffusion evolution equation, and as physical time progresses, use natural diffusion to eliminate noise in the flat areas of the 3D medical volumetric image data, and use normal blocking to maintain the boundary in the edge areas of the 3D medical volumetric image data, until the diffusion evolution equation reaches a quasi-steady state, and output a preprocessed image volume containing an enhanced gradient vector field.

4. The method for segmenting and reconstructing three-dimensional medical image data according to claim 1, characterized in that, In step S2, the Hessian matrix is ​​constructed for the voxels in the preprocessed image volume, and the specific operation for calculating the eigenvalues ​​of the Hessian matrix is ​​as follows: A discrete physical observation aperture sequence containing multiple Gaussian smoothing scale parameters with different values ​​is defined, and a corresponding three-dimensional Gaussian kernel function is constructed for each Gaussian smoothing scale parameter in the discrete physical observation aperture sequence. A three-dimensional convolution operation is performed on the volume of the preprocessed image using a three-dimensional Gaussian kernel function to generate a multi-scale spatial response field. The square of the Gaussian smoothing scale parameter is calculated, and the square of the Gaussian smoothing scale parameter is defined as the scale normalization factor. The multi-scale spatial response field is multiplied and weighted using the scale normalization factor to generate a scale-normalized multi-scale spatial response field. The second-order partial derivative of each voxel in the scale-normalized multi-scale spatial response field in the Cartesian coordinate system is calculated. A 3x3 second-order symmetric matrix describing the curvature of the local energy surface is constructed using second-order partial derivatives. The 3x3 second-order symmetric matrix is ​​defined as the Hessian matrix. Eigenvalue decomposition is performed on the Hessian matrix to obtain three Hessian matrix eigenvalues ​​characterizing the principal curvature. The three Hessian matrix eigenvalues ​​are then sorted in ascending order according to their absolute values. The Hessian matrix eigenvalue with the smallest absolute value among the three Hessian matrix eigenvalues ​​is defined as the first principal curvature along the structural axis, the Hessian matrix eigenvalue with the largest absolute value among the three Hessian matrix eigenvalues ​​is defined as the third principal curvature along the cross-sectional direction, and the Hessian matrix eigenvalue with the middle absolute value among the three Hessian matrix eigenvalues ​​is defined as the second principal curvature.

5. The method for segmenting and reconstructing three-dimensional medical image data according to claim 4, characterized in that, In step S2, the specific operation of removing tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix and extracting regions exhibiting gray-level energy convergence characteristics as the initial seed point set is as follows: Calculate the absolute value of the first principal curvature and define the absolute value of the first principal curvature as the minimum curvature intensity term. The minimum curvature intensity term is used to ensure that the voxel has significant curvature in all directions to eliminate planar structures. Calculate the ratio of the absolute value of the first principal curvature to the absolute value of the third principal curvature. Define the ratio of the absolute value of the first principal curvature to the absolute value of the third principal curvature as the geometric isotropic factor. The geometric isotropic factor is used to quantify the degree to which the local geometric structure approximates a sphere in order to eliminate tubular structures. Define a polarity gate operator, which is used to detect the sign consistency of the three Hessian matrix eigenvalues. When it is necessary to extract bright lesions, configure the polarity gate operator to output a non-zero value when all three Hessian matrix eigenvalues ​​are negative. When it is necessary to extract dark lesions, configure the polarity gate operator to output a non-zero value when all three Hessian matrix eigenvalues ​​are positive. The product of the minimum curvature intensity term, the geometric isotropic factor, and the polarity gate operator is calculated, and the product is defined as the structural saliency response value of the voxel under the Gaussian smoothing scale parameter. The maximum value of the structural significance response of the voxel under all Gaussian smoothing scale parameters is selected, and the maximum value is defined as the multi-scale convergence probability field. The multi-scale convergence probability field represents the probability that the voxel belongs to the core of the lesion. Calculate the non-zero mean and non-zero standard deviation of the multi-scale convergent probability field, calculate the sum of the products of the non-zero mean, non-zero standard deviation and confidence factor, and define the sum as the adaptive truncation threshold. Regions with values ​​greater than the adaptive truncation threshold in the multi-scale convergence probability field are extracted, and a binarized mask is generated. The binarized mask is defined as the initial seed point set.

6. The method for segmenting and reconstructing three-dimensional medical image data according to claim 1, characterized in that, In step S3, an initial level set function is constructed with the initial seed point set as the zero level set region, and a curvature inverse braking term is constructed. The specific operation of constructing the gradient driving term based on the enhanced gradient vector field in step S1 is as follows: Calculate the shortest Euclidean distance from each voxel in 3D space to the edge of the initial seed point set. Construct a symbolic distance field based on the shortest Euclidean distance. The symbolic distance field is defined to take negative values ​​inside the initial seed point set, positive values ​​outside the initial seed point set, and zero values ​​at the edge of the initial seed point set. Define the symbolic distance field as the initial level set function. Define a narrow band region where the initial level set function is close to zero. Update the initial level set function only within the narrow band region. Calculate the vector magnitude of the enhanced gradient vector field output in step S1. Construct an edge indicator scalar function based on the vector magnitude of the enhanced gradient vector field. Set the edge indicator scalar function to approach zero at strong edges where the vector magnitude of the enhanced gradient vector field is large, and set the edge indicator scalar function to approach one at flat regions where the vector magnitude of the enhanced gradient vector field is small. Among them, the vector advection traction term drives the zero-level set interface to adhere to the anatomical structure boundary along the direction of the enhanced gradient vector field, while the scalar expansion thrust term simulates the fluid injection pressure driving the zero-level set interface to expand outward.

7. The method for segmenting and reconstructing three-dimensional medical image data according to claim 6, characterized in that, In step S3, an evolution equation for the level set, consisting of a curvature inverse braking term and a gradient driving term, is established. The initial level set function is iteratively evolved using this equation. Local evolution stops when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function causes the curvature inverse braking term to be greater than the gradient driving term. The specific operation for outputting the evolved level set function is as follows: A partial differential evolution equation following the Hamiltonian Jacobi form is established. The partial differential evolution equation describes the dynamic deformation process of the initial level set function with physical time. A regularized Dirac function is introduced into the partial differential evolution equation, and the evolution force is restricted to the physical location of the zero level set interface by the regularized Dirac function. In the iterative evolution process, for the flat internal region of the anatomical structure, the scalar expansion thrust term dominates, driving the zero-level set interface to grow outward. For the strong edge regions of the anatomical structure, the edge indicator scalar function value approaches zero, causing the evolution rate to approach zero, locking the zero-level set interface on the boundary of the anatomical structure. For the weak boundary adhesion regions of the anatomical structure, the narrow channel causes the local average curvature value to increase sharply, so that the contraction force generated by the curvature reverse braking term exceeds the expansion force generated by the scalar expansion thrust term, causing the local evolution rate to become negative, forcing the zero-level set interface to automatically shrink back and stop growing. The spatial derivative of the initial level set function is discretized using the central difference method, and the temporal derivative of the initial level set function is discretized using the forward Euler method. The time step is dynamically adjusted according to the numerical stability condition, which requires that the time step be less than the minimum value of the ratio of the spatial grid size to the maximum evolution rate. The iterative evolution is continuously executed until the partial differential evolution equation reaches a steady state or the total energy change rate is lower than a preset threshold. The final steady-state scalar field is output as the evolved level set function, and the zero level set interface of the evolved level set function is extracted as the surface of the segmented three-dimensional anatomical structure.

8. The method for segmenting and reconstructing three-dimensional medical image data according to claim 1, characterized in that, In step S4, the evolved level set function output from step S3 is received. The specific operation of extracting the initial mesh vertices and initial vertex connection relationships from the zero level set interface of the evolved level set function is as follows: The zero level set interface in the evolved level set function output by step S3 is parsed using the moving cube isosurface extraction algorithm. The zero level set interface is parsed to obtain the discrete initial mesh vertex set and the initial vertex connection relationship. The initial mesh vertex set contains the initial mesh vertex coordinate information defined in three-dimensional space, and the initial vertex connection relationship contains the triangle mesh list that defines the physical topological adjacency attribute between the initial mesh vertices. The specific operation for calculating the data attachment energy term for the initial grid vertices based on step S1 of enhancing the gradient vector field is as follows: construct a trilinear interpolation operator for the enhanced gradient vector field defined on the discrete voxel grid points output by step S1, and use the trilinear interpolation operator to establish a mapping relationship from the discrete data space to the continuous physical space, so as to transform the discrete enhanced gradient vector field into a continuous gradient potential energy density field. The data-attached energy term physically constructs a continuous gravitational potential well. The energy valley of the gravitational potential well strictly corresponds to the gradient maximum position of the real physical boundary of the anatomical structure. The gravitational potential well generates a gradient adsorption force pointing towards the gradient maximum position. The gradient adsorption force drives the initial mesh vertices to break through the voxel resolution limit. The gradient adsorption force drives the initial mesh vertices to automatically slide towards the position of the local gradient energy maximum in the continuous physical space with sub-voxel precision.

9. The method for segmenting and reconstructing three-dimensional medical image data according to claim 8, characterized in that, In step S4, the elastic deformation energy term for the initial mesh vertices is calculated based on the initial vertex connection relationship. A variational energy functional is constructed, consisting of a weighted sum of the data attachment energy term and the elastic deformation energy term. The optimized spatial coordinates of the initial mesh vertices are calculated by solving for the minimum value of the variational energy functional. The specific operation of generating a 3D mesh model based on the optimized spatial coordinates is as follows: The average curvature vector of the initial mesh vertices is calculated using the discrete Laplace Beltrami operator, and an area-weighted normal vector is introduced to prevent numerical singularities when the initial mesh vertices degenerate. Tangential tension optimizes the shape quality and vertex distribution uniformity of the triangular mesh in the initial vertex connection relationship, while strictly maintaining the volume of the 3D mesh model; Normal pressure eliminates high-frequency noise on the initial mesh vertex surface; A variational energy functional containing data attachment energy and elastic deformation energy is constructed. The variational energy functional is iteratively solved using the semi-implicit Euler method until the maximum displacement of the initial grid vertex is less than the preset thermodynamic convergence threshold. The output is the optimized spatial coordinates of the initial mesh vertex when it reaches mechanical equilibrium. The initial vertex connection relationship remains unchanged. The coordinates of the initial mesh vertex are updated to the optimized spatial coordinates to generate a 3D mesh model with sub-voxel level geometric accuracy and smooth topology.

10. A segmentation and reconstruction system for three-dimensional medical image data, employing the segmentation and reconstruction method for three-dimensional medical image data as described in any one of claims 1-9, characterized in that, include: The system includes an image preprocessing module, a seed point extraction module, a level set evolution module, and a mesh model generation module. The image preprocessing module receives three-dimensional medical volumetric image data, calculates the gradient information of voxels in the three-dimensional medical volumetric image data to construct a structure tensor field, performs eigenvalue decomposition on the structure tensor field to obtain structure tensor eigenvalues, constructs an anisotropic diffusion impedance field based on the structure tensor eigenvalues, substitutes the anisotropic diffusion impedance field into the diffusion equation to perform directional smoothing processing on the three-dimensional medical volumetric image data, and outputs a preprocessed image volume containing an enhanced gradient vector field. The seed point extraction module receives the preprocessed image volume output by the image preprocessing module, constructs a Hessian matrix for the voxels in the preprocessed image volume, calculates the eigenvalues ​​of the Hessian matrix, removes tubular and sheet-like structures based on the eigenvalues ​​of the Hessian matrix, and extracts regions exhibiting gray-scale energy convergence characteristics as the initial seed point set. The level set evolution module receives the initial seed point set output by the seed point extraction module, constructs an initial level set function with the initial seed point set as the zero level set region, constructs a curvature inverse braking term, constructs a gradient driving term based on the gradient vector field enhanced by the image preprocessing module, establishes a level set evolution equation composed of the curvature inverse braking term and the gradient driving term, performs iterative evolution on the initial level set function using the level set evolution equation, stops local evolution when the contraction force generated by the local average curvature at the zero level set interface of the initial level set function makes the curvature inverse braking term greater than the gradient driving term, and outputs the evolved level set function. The mesh model generation module receives the evolved level set function output by the level set evolution module, extracts the initial mesh vertices and initial vertex connections from the zero level set interface of the evolved level set function, calculates the data attachment energy term for the initial mesh vertices based on the enhanced gradient vector field of the image preprocessing module, calculates the elastic deformation energy term for the initial mesh vertices based on the initial vertex connections, constructs a variational energy functional composed of the weighted sum of the data attachment energy term and the elastic deformation energy term, calculates the optimized spatial coordinates of the initial mesh vertices by solving the minimum value of the variational energy functional, and generates a three-dimensional mesh model based on the optimized spatial coordinates.

Images