Heart three-dimensional polygon mesh construction method and system

Abstract

The invention relates to the technical field of medical image calculation, in particular to a heart three-dimensional polygon mesh construction method and system, and the method comprises the steps: firstly obtaining and registering an image, generating an inductor density and a signed distance field based on boundary features, and fitting anatomic anchor points; determining a patient specific three-dimensional domain through mass conservation transmission and homeomorphic deformation; constructing pull-back measurement in the domain, solving three groups of harmonic functions, extracting contour surfaces according to common area constraint, and intersecting to form a network; deformation and harmonic functions are taken as joint unknown quantities, synchronous optimization is carried out in combination with pixel consistency, physical consistency and topological consistency, and the density of the inducer is re-weighted by layer area deviation feedback until convergence. According to the method, the surface net with consistent anatomy, controlled topology and balanced layer is output, and the method is adaptive to diagnosis, treatment and simulation application.

Description

Technical Field

[0001] This invention relates to the field of medical image computing technology, and in particular to a method and system for constructing a three-dimensional polygonal mesh of the heart. Background Technology

[0002] Image-driven personalized cardiovascular modeling has become the fundamental data carrier for diagnostic quantification, device selection, preoperative planning, and numerical simulation. Mesh quality directly affects multi-center comparability and clinical implementation efficiency. Existing step-by-step workflows rely heavily on intensity registration and direct isosurface reconstruction, which are prone to volume drift and topological folding, neglect deformation-induced anisotropy, have insufficient sampling in key areas, and lack global feedback between stages, resulting in unstable quality and high costs for manual remeshing. Summary of the Invention

[0003] To address the numerous problems existing in the prior art, this invention provides a method and system for constructing a three-dimensional polygonal mesh for the heart. This invention uses mass conservation to transfer and align the template volume density and the induced volume density, and obtains homeomorphic deformation; it constructs a pull-back metric in the patient-specific three-dimensional domain, solves three sets of harmonic functions consistent with the anatomical direction, and extracts isosurfaces to form a mesh based on co-area constraints; it uses pixel consistency, physical consistency, and topological consistency to form a joint target closed-loop optimization, thereby obtaining a high-quality surface mesh that is anatomically aligned, topologically controlled, and layer-balanced.

[0004] This specification provides one or more embodiments of a method for constructing a three-dimensional polygonal mesh of the heart, including the following steps: Acquire and register medical images, construct inducible body density and signed distance field based on boundary features, and fit the set of anatomical anchor points; Mass conservation transport is calculated between template body density and inducer body density to obtain transport mapping, homeomorphic deformation is calculated and constrained at the anatomical anchor set, and the patient-specific cardiac three-dimensional domain is determined by combining pixel consistency and physical consistency. A pullback metric is constructed in the three-dimensional domain of the patient's specific heart. Three sets of harmonic function fields are solved and aligned with the anatomical direction. Isosurfaces are extracted based on co-area constraints and their intersections are used to generate a polygonal surface mesh. Using homeomorphic deformation and three sets of harmonic function fields as joint unknowns, a target including pixel consistency, physical consistency and topological consistency is constructed. The polygonal surface mesh is updated and reconstructed synchronously. The density of the inducer is adjusted based on the layer area deviation feedback until convergence, and the three-dimensional polygonal surface mesh of the heart is output.

[0005] According to one or more embodiments of this specification, acquiring and registering medical images includes performing spatial registration and intensity normalization, extracting boundary features based on gradient structure tensor and curvature estimation, constructing inducible body density from the boundary features through kernel smoothing, calculating the signed distance field from the topologically purified initial boundary using the fast travel method, and the anatomical anchor set includes the valve annulus closure curve, the apex of the heart, the opening of the great vessels, and the master curve of the interventricular septum.

[0006] According to one or more embodiments of this specification, calculating the mass-conserving transport between the template volume density and the inducer volume density includes solving for the optimal transport using entropy regularization, with the cost function being the squared Euclidean distance.

[0007] According to one or more embodiments of this specification, the method for calculating homeomorphic deformation includes calculating homeomorphic deformation based on an exponential mapping of a static velocity field, applying positional constraints at the set of anatomical anchor points, setting non-crossable boundary constraints in the atrioventricular septal region, and setting positive constraints on the Jacobian determinant to maintain the mapping as homeomorphic.

[0008] According to one or more embodiments of the method described in this specification, pixel consistency is calculated by the reprojection error of differentiable rendering, and physical consistency is calculated by the residuals of the forward solution and the adjoint solution of the wave equation. The gradients of the two are backpropagated along the mapping relationship of homeomorphism to update the homeomorphism and determine the patient-specific cardiac three-dimensional domain.

[0009] According to one or more embodiments of this specification, constructing a pull-back metric includes pulling back an Euclidean metric through homeomorphic deformation to obtain a metric tensor, applying Dirichlet boundary conditions to three sets of harmonic function fields at the apex and base, applying Dirichlet boundary conditions to the interventricular septum and applying Neumann boundary conditions to the lateral wall, applying Dirichlet boundary conditions to the endocardium and epicardium, and applying Neumann boundary conditions to the opening loop of the great vessels.

[0010] According to one or more embodiments of this specification, the method for extracting isosurfaces based on co-area constraints includes constructing weights based on curvature and anatomical region indicator functions, increasing weights in the valve ring and outflow channel regions to form local densification according to the weight allocation level, and generating polygonal surface mesh units aligned with the anatomical direction on the intersection line of isosurfaces.

[0011] According to one or more embodiments of this specification, the synchronous update includes constructing a block linear system with Karush-Kuhn-Tucker constraints as joint unknowns using homeomorphic deformation and three sets of harmonic function fields. The constraints include pixel consistency constraints, physical consistency constraints and topological consistency constraints. An iterative linear solver is used to solve the block linear system to achieve synchronous update and reconstruction of the polygonal surface mesh.

[0012] According to one or more embodiments of this specification, adjusting the induced body density based on layer area deviation feedback includes converting the layer area deviation into a density reweighting factor, reweighting the induced body density and using it as the target density in the next mass conservation transfer, and repeating synchronous updates and reconstructions until the polygonal surface mesh meets the preset mesh quality requirements and topology consistency requirements.

[0013] Some embodiments of this specification also provide a cardiac three-dimensional polygonal mesh construction system for implementing the described cardiac three-dimensional polygonal mesh construction method, the system comprising: The image modeling module is used to acquire and register medical images, construct inducible body density and signed distance fields based on boundary features, and fit the set of anatomical anchor points. The transport deformation module is used to calculate the mass conservation transport between the template body density and the inducer body density to obtain the transport mapping, calculate the homeomorphic deformation and constrain it at the anatomical anchor point set, and combine pixel consistency and physical consistency to determine the patient-specific cardiac three-dimensional domain. The metric meshing module is used to construct pullback metrics in the three-dimensional domain of the patient's specific heart, solve three sets of harmonic function fields and align them with the anatomical direction, extract isosurfaces based on co-area constraints and generate polygonal surface meshes by their intersection; The closed-loop optimization module is used to construct a target that includes pixel consistency, physical consistency and topological consistency using homeomorphic deformation and three sets of harmonic function fields as joint unknowns. It synchronously updates and reconstructs the polygonal surface mesh, adjusts the inducing body density based on layer area deviation feedback until convergence, and outputs a three-dimensional polygonal surface mesh of the heart.

[0014] Compared with the prior art, the advantages and beneficial effects of the present invention are as follows: By using mass conservation transport and homeomorphic deformation coupling constraints, stable alignment and topology preservation were achieved under conditions of strength difference and large deformation.

[0015] By aligning the pullback metric with the orientation of the three sets of harmonic fields, parameterized stratification consistent with the anatomical axis and trans-wall anisotropy fidelity were achieved.

[0016] By combining common area constraints with anatomical region weights, we achieved balanced layer area and adaptive densification of key areas such as valves and outflow channels.

[0017] By simultaneously optimizing the joint objectives of pixel consistency, physical consistency, and topological consistency, we achieved coordinated convergence of image alignment, prior physical rules, and topological rules.

[0018] By using the induced volume density reweighting with layer area deviation as feedback and outer layer closed-loop iteration, cross-stage consistency maintenance and controllable convergence of mesh quality are achieved. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the execution flow of the method of the present invention; Figure 2 This is a structural block diagram of the system of the present invention. Detailed Implementation

[0020] The embodiments of the present disclosure will now be described with reference to the accompanying drawings. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the disclosure. In the following detailed description, numerous specific details are set forth to provide a thorough understanding of the embodiments of the present disclosure for ease of explanation. However, it will be apparent that one or more embodiments may be practiced without these specific details. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concepts of the present disclosure.

[0021] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit this disclosure. The terms “comprising,” “including,” etc., as used herein indicate the presence of the stated features, steps, operations, and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, or components.

[0022] All terms used herein (including technical and scientific terms) have the meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein are to be interpreted in a manner consistent with the context of this specification, and not in an idealized or overly rigid way.

[0023] In this invention, polygon mesh construction is considered a "geometric compilation" process that transforms voxel-level image information into computationally usable geometric representations. Its goal is to describe continuous cardiac structures with a finite and consistent set of facets, while simultaneously satisfying four constraints: geometric fidelity, topological correctness, parametric field usability, and cross-step consistency. Specific requirements include: maintaining shape and position consistency at anatomical anchor points and key curves; maintaining slice order and monotonicity in transmural directions; achieving controlled densification and facet quality constraints in the vicinity of the interventricular septum, outflow tract, and valve annulus; and ensuring the repeatability of normal, curvature, and area statistics in subsequent simulations and measurements. Traditional workflows first extract isosurfaces directly from images and then re-mesh them, which struggles to simultaneously satisfy the above constraints and easily leads to slice imbalance, folding, and cross-step inconsistencies. This invention revolves around the core concept of polygonal mesh construction. It proposes a patient-specific domain determination mechanism based on mass conservation and homeomorphism. Combining pullback metric and directional parameter fields of three sets of harmonic functions, it generates isosurfaces under co-area constraints and they intersect to form a mesh. Finally, it updates the deformation and parameter fields synchronously with a closed-loop target consisting of pixel consistency, physical consistency and topological consistency, and outputs an anatomically aligned, topologically controlled, and layer-equilibrium polygonal surface mesh that can be directly used for calculation and visualization.

[0024] like Figure 1As shown, a method for constructing a three-dimensional polygonal mesh for the heart includes the following steps: Acquire and register medical images, construct inducible body density and signed distance field based on boundary features, and fit the set of anatomical anchor points; Acquiring and registering medical images includes spatial registration and intensity normalization, extracting boundary features based on gradient structure tensor and curvature estimation, constructing inducible body density from boundary features through kernel smoothing, calculating the signed distance field from the topologically purified initial boundary using the fast travel method, and the set of anatomical anchor points includes the valve annulus closure curve, the apex of the heart, the opening loop of the great vessels, and the master curve of the interventricular septum.

[0025] This embodiment focuses on the preprocessing and geometric constraint data construction of cardiac images. The goal is to obtain the inducible body density, signed distance field, and anatomical anchor set under a unified coordinate and intensity system, which can serve as stable inputs for subsequent deformation estimation and network formation.

[0026] First, import cardiac computed tomography (CT) or magnetic resonance imaging (MRI) images, unifying them to fixed voxel spacing and world coordinates. Spatial registration is then performed, initially using rigid or affine registration based on mutual information to align the patient image to the reference anatomy. Next, restricted non-rigid registration is performed under cardiac coarse segmentation mask guidance to correct local deformations while avoiding overstretching. After registration, intensity normalization is performed: bias field correction is applied to the MRI images to remove intensity inhomogeneities; intensity recalibration or histogram matching is performed within the statistical range of myocardium and blood cavities to suppress cross-device differences. Subsequently, three-dimensional edge-preserving denoising is performed, preferentially using anisotropic diffusion or nonlocal means to preserve boundary gradients.

[0027] The inducible body density and signed distance field are constructed based on boundary features. To robustly extract boundary candidates, the gradient structure tensor is first calculated to obtain the dominant normal direction, and the boundary confidence is evaluated accordingly. Simultaneously, high curvature regions are identified using curvature estimation. Topological purification of the initial boundary is performed using connected component analysis, void filling, and small branch removal to ensure that the connectivity between the cardiac chambers and myocardium conforms to anatomical common sense. Using the purified binary boundary as a seed, the signed distance field is solved using a fast travel method, defining the cardiac chamber direction as negative and the external myocardial direction as positive, with zero-faces corresponding to geometric interfaces such as the endocardium and endocardium. To obtain a continuously optimizable target distribution within the volume domain, high-confidence boundary points are kernel-smoothed diffused within the volume, with the bandwidth adaptively adjusted based on local curvature and point spacing. Boundary confidence is used as a weight to obtain the inducible body density. This density provides sufficient contrast in the boundary neighborhood and gradually decays within and outside the cavity and myocardium, facilitating numerical stability of downstream mass-conserving transport.

[0028] The fitting of the anatomical anchor set focuses on four elements: the valve annulus closure curve, the apex, the great vessel opening loops, and the interventricular septum master curve. The extraction of the valve annulus closure curve employs a slice sequence method: continuous orthogonal slices are made near the estimated valve normal direction, strong annular boundaries are detected, and robust splines are fitted. These slices are then concatenated to form a closed curve. To avoid annulus breakage, a regularization term from the distance field is introduced into low-contrast slices to maintain curve smoothness. The apex is determined based on a combined criterion of centerline and wall thickness: the left ventricular centerline is obtained using the distance field framework or constrained shortest path, and the apex is located far from the base where the local wall thickness tends to be minimal. The extraction of the great vessel opening loops is performed within the outflow tract voxel zone. Near-circular or elliptical boundaries are detected through orthogonal slices, and circles or ellipses are fitted and projected back into three-dimensional space to form the opening loops. If necessary, the directional continuity of adjacent slices is used as a constraint. The extraction of the interventricular septum master curve involves curve tracing within the minimum thickness zone, and lateral wall suppression terms are used to avoid deviation to the free wall, ultimately outputting a continuous curve consistent with the anatomy.

[0029] To ensure feasibility, all results are output in a unified coordinate and data format: the inducible body density is presented as a three-dimensional scalar field; the signed distance field is also presented as a three-dimensional scalar field with a zero-level set grid; the anatomical anchor set is recorded as a sequence of curve control points and curve types, along with compatibility verification indices between the curves and the distance field. Quality control recommendations include: checking the closure and single connectivity of the valve annulus and ventral opening loops, the consistency between the apex position and the end of the central axis, the overlap between the zero-level set of the signed distance field and the boundary confidence region, and the monotonicity of the inducible body density on both sides of the boundary. Through this process, the basic geometric and statistical quantities required for subsequent deformation estimation and polygon networking can be stably obtained without complex theoretical derivations, facilitating consistent access and reproduction of images from different institutions and multiple imaging modalities.

[0030] Mass conservation transport is calculated between template body density and inducer body density to obtain transport mapping, homeomorphic deformation is calculated and constrained at the anatomical anchor set, and the patient-specific cardiac three-dimensional domain is determined by combining pixel consistency and physical consistency. This embodiment illustrates the mass conservation and transfer between template volume density and induced volume density, the constrained solution of homeomorphic deformation, and the joint criterion of pixel consistency and physical consistency, used to determine the patient-specific cardiac three-dimensional domain. The workflow is designed for engineering implementation, employing a multi-scale and stepwise convergence strategy to avoid high sensitivity to image quality and initial values.

[0031] First, the template volume density and induced volume density are uniformly represented on a voxel mesh. Mass conservation transfer is performed using a pyramid resolution, from coarse to fine. The coarse scale obtains the global correspondence, while the fine scale supplements local details. The initial values ​​of the transfer mapping can be given by centroid matching and histogram alignment. Subsequently, a stable mapping is obtained through optimal transfer iteration with entropy regularization. To maintain consistency with subsequent geometry, the transfer mapping only provides large-scale mass allocation, while homeomorphic deformation is responsible for geometric continuity and local invertibility.

[0032] The core optimization objective adopts a single joint formula, and the joint formula calculation expression is as follows: ; in, ; The density of the template volume. For the density of the inducible body, For transport mapping, It is a homeomorphic transformation. For static velocity field, For the set of anatomical anchor points, This is the reprojected image after deformation. To observe the images, For physical consistency residuals, and and These are the weighting coefficients.

[0033] In the specific implementation, an initial solution for the transport mapping is first obtained at a coarse scale by updating the iterative ratio of optimal transport using entropy regularization, and neighborhood pruning is performed on abnormally large displacement regions to ensure the mapping is monotonic. Then, homeomorphic deformation is calculated, and an exponential mapping is achieved using the scalar integral and scale-squaring method of the static velocity field to ensure continuous and reversible deformation. To implement anatomical constraints, key point positions are bound by hard constraints at the anatomical anchor point set, and an impassable zone is set near the atrioventricular septum and valve annulus, combined with a positive determinant penalty at the voxel level to suppress folding. Pixel consistency uses the absolute error of differentiable reprojection, while physical consistency uses the low-frequency response model of myocardial tissue or the residual of a validated propagation model; the gradients of both are backpropagated to the velocity field and transport mapping during optimization.

[0034] The solution strategy employs alternating minimization and trust region step size control: the velocity field is updated with a fixed transport map, and then the transport map is refined with the velocity field fixed again; at the end of each round, determinant and boundary crossing checks are performed, and if constraints are violated, the process is rolled back and the step size is reduced. To improve robustness, anisotropic smoothing is introduced in the neighborhood of the ejection duct and valve annulus, allowing displacement updates only in the tangential direction, while normal updates are restricted. After joint convergence, the template cardiac 3D domain is mapped to the patient space through homeomorphic deformation to obtain the patient-specific cardiac 3D domain; simultaneously, the minimum determinant, anchor point reprojection error, pixel error, and physical residual curves are output as quality control records.

[0035] Through the above steps, under the premise of satisfying mass conservation and topological invariance, the imaging and physiological constraints can be aligned simultaneously to form a stable patient-specific cardiac three-dimensional domain that can be used for subsequent metric construction and networking.

[0036] The calculation of mass-conserving transport between the template volume density and the induced volume density involves solving for the optimal transport using entropy regularization, with the cost function being the squared Euclidean distance.

[0037] This embodiment focuses on the mass conservation and transfer between the template volume density and the induced volume density, producing a transfer map and displacement field as input for subsequent homeomorphic deformation solutions.

[0038] First, the template volume density and the induced volume density are represented on the same 3D voxel mesh, with the voxels having the same size and confined within the effective mask of the heart. The ratio of their total mass is calculated and scaled proportionally to make their total masses equal; a lower limit is set for the minimum voxel values ​​to avoid irreversibility caused by values ​​of zero. A multi-scale pyramid is constructed, and the solution is obtained by solving stepwise from coarse to fine to obtain a stable solution.

[0039] The cost is calculated using the squared Euclidean distance and locally sparsified, retaining only pairs within a certain neighborhood radius to reduce computational cost. Optimal transport is obtained using entropy regularization and solved through iterative scaling. The transport solution expression is given below:

[0040]

[0041]

[0042] in, For template mesh number Individual center coordinates, Patient Grid Individual center coordinates, The square of the Euclidean distance. Here is the entropy regularization coefficient. It is a Gibbs core. Let be the mass vector of the template volume density across each voxel. Let the induced volume density be the mass vector of each voxel. As a proportional vector, For discrete transportation planning, For the first The centroid corresponding point of each template position For displacement vectors, To start from the template voxel center To the patient voxel center The allocated transport mass represents the mass of the transport in a mass-conserving transport system, determined by position. Transport to location The quality share.

[0043] The implementation details are as follows: Initialize with a large regularization coefficient at the coarsest scale, and obtain the initial value by updating the scale 10 to 20 times. When refining step by step, interpolate the values ​​of the previous level proportionally. As initial values, the regularization coefficient is gradually reduced to restore details. After each update, the marginal constraint residual is calculated and used as a convergence criterion; if the residual does not meet the criterion, iteration continues. To avoid numerical underflow, safe division and exponentiation operations are implemented in the logarithmic field for the kernel matrix; the kernel matrix is ​​based on sparse neighborhood storage, combined with batch multiplication to improve efficiency. Pairing of the extracardiac region is directly masked outside the mask to ensure that the transport plan allocates quality only within a reasonable anatomical range.

[0044] After obtaining the transport plan, the centroid corresponding point and displacement vector are calculated according to the transport solution expression. The displacement field is anisotropically smoothed within the myocardial mask, allowing fine-tuning along the wall tangential direction, while the normal direction is constrained with a small step size to suppress boundary crossing. For voxels exhibiting isolated large displacements, neighborhood median replacement is performed, and morphological constraints are applied around the valve annulus and outflow tract to ensure displacement continuity. The displacement field is superimposed onto the template coordinates to form the transport map; simultaneously, the mass conservation error, edge residuals, and maximum displacement amplitude are output as quality control indicators. If the mass conservation error or residual exceeds the threshold, the process returns to the previous scale, increases the number of iterations, or slightly increases the regularization coefficient before retrying.

[0045] Through the above steps, under the premise of ensuring mass conservation and numerical stability, we can obtain the transport mapping and displacement field consistent with the anatomical region, which provides a good initial value for the constrained optimization of homeomorphic deformation and significantly reduces the search difficulty and convergence pressure of subsequent geometric solutions.

[0046] The calculation of homeomorphic deformation includes calculating homeomorphic deformation based on the exponential mapping of the static velocity field, applying positional constraints at the set of anatomical anchor points, setting non-crossable boundary constraints in the atrioventricular septal region, and setting positive constraints on the Jacobian determinant to maintain the mapping as homeomorphic.

[0047] This embodiment provides a solution process for homeomorphic deformation. The inputs are the initial displacement values ​​obtained from the transport mapping, the set of anatomical anchor points, and the inaccessible compartmentalization region determined based on the distance field. The output is a three-dimensional deformation that satisfies connectivity and reversibility and passes the quality and geometric consistency checks.

[0048] The static velocity field is used to generate a continuously reversible deformation flow. The calculation expression is as follows: In three-dimensional coordinates, For a moment Deformation flow, For static velocity field, The target is homeomorphic deformation. In implementation, the scale-flat method is used for numerical exponential mapping, first approximating at a coarse scale, and then refining step by step to ensure numerical stability.

[0049] The velocity field is initialized from the displacement field of the transport map, and initial values ​​are obtained through anisotropic smoothing. Details are preserved along the myocardial tangent, and oscillations are suppressed along the normal. Then, iterative optimization is performed: first, the gradient information of the objective function is formed based on the pixel reprojection error and the physical consistency residual; second, the excessive changes in the velocity field are suppressed based on the regularization term; and third, the velocity field and deformation flow are corrected according to constraints so that the constraints are satisfied after each iteration.

[0050] The anatomical anchor points are implemented using a combination of hard constraints and small-range elastic bands. First, the anchor point correspondence is bound to the deformation end position. After each update of the velocity field, the reprojection error of the anchor point is calculated and the velocity field is adjusted in the local neighborhood to ensure that the anchor point is accurately positioned, while limiting excessive energy diffusion to distant locations to avoid introducing non-physical anatomical deformation.

[0051] The atrioventricular septum is made non-crossable by constructing a normal band using the distance field. Within this band, the unit normal is calculated, and the velocity field is set to zero in the normal direction, retaining only the tangential component. For generated deformable flows, if a crossing indication occurs, the corresponding node is repositioned to the nearest legal side, and normal suppression is strengthened in the next velocity update until the crossing count is zero. This process ensures that the geometric separation between the interventricular septum and the annular region is not disrupted.

[0052] To ensure that the mapping is invertible everywhere, a positive constraint is imposed on the Jacobian determinant, and numerical barriering is applied: For Jacobi determinant, The threshold is positive. The implementation method is as follows: Calculate using finite differences on a voxel mesh. When the local minimum value approaches the threshold, the step size is automatically reduced and the regularization weight is increased; if individual positions are below the threshold, the last update is rolled back and the volume compression component of the velocity field is restricted in that neighborhood until it recovers to above the threshold.

[0053] Numerical integration employs a multi-scale strategy: coarse-scale integration uses larger step sizes and stronger regularization to quickly obtain a globally consistent deformation trend; fine-scale integration gradually reduces the step size to recover details. After each scale is completed, three quality controls are performed: maximum anchor point deviation, count of crossings of uncrossable regions, and the global minimum of the Jacobian determinant. Only when all three indicators simultaneously meet the criteria is the process moved to a finer scale.

[0054] To improve robustness, directional smoothing is applied around the valve annulus and outflow tract, allowing velocity updates only along the annulus and flow directions. Sparse constraints are added near the apex to prevent localized overstretching. The final deformation is a homeomorphic mapping, and the output includes the deformation field, velocity field, anchor point deviation statistics, minimum Jacobian determinant, and a non-crossing verification report. This embodiment can stably generate patient-specific three-dimensional domains while ensuring no crossing of critical structures and no folding, providing a reliable geometric basis for subsequent metric construction and meshing.

[0055] Pixel consistency is calculated from the reprojection error of differentiable rendering, while physical consistency is calculated from the residuals of the forward and adjoint solutions of the wave equation. The gradients of the two are backpropagated along the mapping relationship of homeomorphism to update the homeomorphism and determine the patient-specific cardiac three-dimensional domain.

[0056] This embodiment presents an implementable process for the joint constraints of pixel consistency and physical consistency. It takes the current homeomorphic deformation as input and outputs a patient-specific cardiac 3D domain that satisfies both imaging alignment and physical constraints. Pixel consistency uses a differentiable renderer to reproject the deformed volume data under the original acquisition geometry. Errors are only statistically analyzed within the foreground mask. The rendering kernel uses trilinear interpolation and a fixed volume transfer function setting to ensure differentiability. Physical consistency uses a validated low-frequency wave propagation model as a priori. A simulated signal is constructed through forward simulation and compared with the observed signal. The gradient source is obtained from the adjoint solution through time inversion. The calculation expression for the joint objective is: For differentiable rendering operators, For homeomorphic deformation mapping, To observe the images, The simulated signal is calculated within the deformation domain. For observing signals, For the cardiac computation domain, For signal duration, and The weighting coefficients are used. The gradient of the above objective is passed to the voxel coordinates via the rendering derivative and the physical adjoint field, and then accumulated to the velocity field update via the chain relationship from coordinates to deformation. Homeomorphism constraints and non-crossing constraints remain valid through the projection step.

[0057] The engineering implementation employs pyramid resolution and alternating updates. The first step involves fixing the homeomorphic deformation, calculating pixel reprojection, generating a residual map, and accumulating it by voxel to obtain a pixel-consistent gradient field. The second step involves fixing the pixel terms, running a forward wave propagation simulation to generate a simulated signal, performing time inversion to obtain the adjoint field, and backpropagating the physical residual to the coordinate system to obtain a physically consistent gradient field. The third step involves summing the two types of gradients by weight, clearing the normal component to zero within the inaccessible band, retaining only the tangential update, and then anisotropically smoothing the velocity field. The homeomorphic deformation is then updated using a line search or trust region strategy. At the end of each iteration, the minimum values ​​of anchor point deviation, traverse count, and Jacobian determinant are checked. If any of these indicators fail to meet the criteria, the iteration is backtracked and the step size is reduced.

[0058] To reduce sensitivity to noise, the pixel term assigns higher weights to regions with strong gradients and zero weights to the background; the physics term samples within the myocardium and reduces weights within large blood vessel cavities. The camera geometry of the differentiable renderer is read from the medical image header file, and temporal and spatial sampling use a consistent step size to avoid artifacts. The convergence criterion is that the target descent rate, maximum anchor point deviation, and Jacobian determinant threshold must all be satisfied simultaneously. The final output is a patient-specific cardiac 3D domain along with corresponding error statistics and quality control reports, which are used in subsequent metric construction and network generation steps.

[0059] A pullback metric is constructed in the three-dimensional domain of the patient's specific heart. Three sets of harmonic function fields are solved and aligned with the anatomical direction. Isosurfaces are extracted based on co-area constraints and their intersections are used to generate a polygonal surface mesh. In this embodiment, a pullback metric is constructed within the three-dimensional domain of the patient's specific heart. Three sets of harmonic function fields are solved and aligned with the anatomical direction. Then, isosurfaces are extracted under co-area constraints and their intersections generate a polygonal surface mesh, which serves as the basis for subsequent mesh optimization and calculation.

[0060] The input is a patient-specific cardiac 3D domain and homeomorphic deformation mapping. First, the deformation Jacobian and gradient are calculated on a voxel or tetrahedral mesh to construct a pullback metric tensor. The expression for calculating the pullback metric tensor is: ; To pull back the metric tensor, For homeomorphic deformation mapping, This is the spatial gradient operator. During implementation... Pruning is performed based on eigenvalue decomposition, limiting the minimum and maximum eigenvalues ​​to avoid numerical ill-conditioning, and edge-preserving smoothing is performed in the neighborhood of the elements to ensure that anisotropy along the myocardial wall thickness direction is preserved while the tangential direction is not excessively blurred.

[0061] The patient-specific cardiac 3D domain was discretized into a quality-controlled tetrahedral mesh. Inferior elements (too small in volume or too large in aspect ratio) were removed, and the mesh was appropriately densified near the valve annulus and the opening of the great vessels to improve geometric fidelity. The pullback metric tensor was assigned to the mesh using a center-averaged method and used as the coefficient field for subsequent equations.

[0062] The three sets of harmonic function fields are aligned with three principal anatomical directions: along the apex to the base, along the interventricular septum to the lateral wall, and along the endocardium to the epicardium. The harmonic functions are expressed using a harmonic equation under the pullback metric, with the following expression: The same formula applies to the other two sets of harmonic function fields. This represents the harmonic function field. Boundary conditions are set according to anatomical constraints: constant boundaries are set at the apex and base to determine the scaling in the transbasal direction; complementary constant or flux boundaries are set at the interventricular septum and lateral walls to fix the transverse orientation; complementary boundaries are set at the endocardium and epicardium to ensure transmural monotonicity. For the valve annulus closure curve and the great vessel opening loop, the harmonic function values ​​are smoothly transitioned along the circumferential direction by curve interpolation to avoid discontinuities in the junctional region. The linear equation system is solved using a conjugate gradient algebraic multigrid with preconditions, and the convergence criterion is that the residual descent and boundary error simultaneously satisfy a threshold.

[0063] After completing the three sets of harmonic function fields, alignment and orthogonality correction are performed: the three-directional gradient is calculated within each cell; if the included angle is too small or reversed, least squares orthogonalization is used for local correction, and the result is written back to the function value to maintain global continuity. Subsequently, layer allocation is generated based on the common area constraint: using the trans-wall harmonic function as the main control variable, several thresholds are uniformly sampled within the interval, and the initial layers are obtained by extracting isosurfaces. After calculating the area of ​​each layer, the thresholds are slightly adjusted through one-dimensional monotonic mapping to make the area of ​​each layer close to the target allocation; the threshold update adopts binary search or Newton-style stepping to ensure that the layer order remains unchanged and the monotonicity is not violated.

[0064] Isosurface extraction employed a triangular polygonization algorithm (isosurface cube method for voxel data, triangulation method for tetrahedral data), yielding three families of isosurfaces. The three families of isosurfaces were intersected pairwise to obtain grid lines, and their intersections formed corner points. Polygonal surface mesh elements were generated by the closed segmentation of the grid lines on the facets. To obtain usable polygonal surface meshes, small-angle facets were locally re-divided, short sides were collapsed, and serrated boundaries were smoothed to preserve their shape. Topological operations were prohibited in the neighborhood of the valve annulus, the main line of the interventricular septum, and the opening loop; only tangential smoothing was allowed to maintain the position and shape of the anatomical curves.

[0065] Quality control includes: checking the monotonicity of the transwall harmonic function from the inner membrane to the outer membrane, checking the target area deviation of each layer, and checking the manifold and porosity of the polygonal surface mesh; triggering automatic backtracking for unsatisfactory items, fine-tuning the threshold, or re-extracting the mesh locally. The final output is a polygonal surface mesh that is consistent with the anatomical direction, has controlled layer area, and is topologically correct, along with the corresponding layer threshold and quality report, which can be directly used for subsequent closed-loop optimization and mesh application.

[0066] The construction of the pull-back metric involves pulling back the Euclidean metric through homeomorphic deformation to obtain the metric tensor. Three sets of harmonic function fields are subjected to Dirichlet boundary conditions at the apex and base, Dirichlet boundary conditions at the interventricular septum and Neumann boundary conditions at the lateral wall, Dirichlet boundary conditions at the endocardium and epicardium, and Neumann boundary conditions at the opening of the great vessels.

[0067] Within the patient-specific cardiac three-dimensional domain, the pullback metric tensor is first calculated based on the aforementioned allomorphic deformation mapping, with the calculation relationship described in the pullback metric tensor calculation expression above. To ensure numerical stability, eigenvalue pruning is performed on the pullback metric tensor, and edge-preserving smoothing is applied to adjacent cells, thereby suppressing pathological conditions while maintaining anisotropic information in the transmural direction.

[0068] The domain is discretized into a quality-controlled tetrahedral mesh, eliminating elements with excessively small volumes or excessively high aspect ratios. Appropriate densification is applied to the neighborhoods of the valve annulus and the openings of major blood vessels to improve geometric fidelity. The pullback metric tensor is written into the discretization equations as a field of element coefficients. The three sets of harmonic function fields correspond to the three principal anatomical directions: apex to basement, interventricular septum to lateral wall, and endocardium to epicardium, respectively. Their governing equations are as follows: The boundary conditions were set as follows: Dirichlet boundary conditions were applied at the apex and basal region to determine the transbasal scaling; Dirichlet boundary conditions were applied at the interventricular septum and Neumann boundary conditions were applied at the lateral walls to limit lateral orientation and flux; Dirichlet boundary conditions were applied at the endocardium and epicardium to ensure transwall monotonicity; and Neumann boundary conditions were applied around the openings of large blood vessels to avoid non-physical flux. The linear system was solved using conjugate gradients combined with algebraic multigrid preconditioning, and the convergence criterion was that both residual descent and boundary error simultaneously met a threshold.

[0069] After the solution is obtained, direction alignment and orthogonality correction are performed: the gradient directions of the three sets of harmonic fields are calculated for each element. If local anomalies such as excessively small included angles or reversed directions occur, least squares orthogonalization is used for local correction, and the results are written back to the node values ​​to maintain global continuity and orthogonality. Subsequently, layer allocation is performed based on the common area constraint: the trans-wall harmonic field is used as the main control variable. First, the threshold is uniformly sampled within the value interval to extract the initial isosurface and count the area of ​​each layer. Then, the threshold is slightly adjusted through one-dimensional monotonic mapping to make the layer area close to the desired allocation. The threshold update maintains the layer order and monotonicity without being destroyed.

[0070] Isosurface extraction employs a polygonization algorithm matched to the mesh type, yielding three families of isosurfaces. Pairwise intersections of these isosurfaces generate mesh lines, and the intersections of the three families form corner points, thus completing in-plane segmentation to obtain polygonal surface mesh elements. Local re-meshing and shape-preserving smoothing are performed on small-angle patches and serrated boundaries, while edge collapse is applied to short sides. Topological modifications are prohibited in the neighborhood of the valve annulus, the master curve of the interventricular septum, and the opening loop; only tangential smoothing is allowed to ensure that the position and shape of key anatomical curves are not destroyed.

[0071] Quality control includes: checking the monotonicity of the transmural harmonic field in the direction from the endocardium to the epicardium; checking that the horizontal sets in the direction from the interventricular septum to the lateral wall do not cross; and checking that the horizontal sets in the direction from the apex to the base do not fold back. Local re-solution or boundary re-labeling is triggered for any non-satisfied terms. The final output is a polygonal surface mesh that is consistent with the anatomical direction, has controlled slice area, and is topologically correct, along with the corresponding slice threshold and quality report.

[0072] The extraction of isosurfaces based on co-area constraints includes constructing weights based on curvature and anatomical region indicator functions. According to the weight allocation level, the weights are increased in the valve ring and outflow channel regions to form local densification, and polygonal surface mesh units aligned with the anatomical direction are generated on the intersection line of the isosurfaces.

[0073] This embodiment focuses on the layer extraction and meshing stage. The goal is to generate a family of isosurfaces aligned with the anatomical direction within the patient-specific three-dimensional cardiac domain based on co-area constraints, and to locally densify the isosurfaces in the valve annulus and large vessel outflow tract regions. Finally, polygonal surface mesh units are constructed from the intersection lines of the isosurfaces.

[0074] First, a hierarchical weighting system is constructed. Multi-scale second-order differences are performed on the transwall harmonic function field to obtain voxel-level curvature indices, and edge-preserving smoothing is applied near the anatomical boundaries. A comprehensive weighting field is then formed by combining the regional indicator masks of the valve annulus and the outflow tract of the large vessels. The expression for the comprehensive weighting field is as follows: For curvature normalization, This is an indicator function for the anatomical region. and These are weighting coefficients. The aforementioned weighting fields take larger values ​​in high curvature and critical regions to guide subsequent layer density allocation.

[0075] Based on the concept of common area, the distribution is performed in the trans-wall direction. Let the trans-wall harmonic function be denoted as a scalar field, and the weighted cumulative volume function is calculated as follows: This represents the weighted cumulative volume from the intima side to the threshold surface. For transwall harmonic functions, For the threshold, As weight, Let be the spatial volume element. Given the number of layers, distribute the total cumulative volume evenly and solve for the condition that... A set of monotonic thresholds The solution employs a monotonic bisection or secant method, numerically checking the threshold order and monotonicity after each iteration. In the valve annulus and large vessel outflow tract regions, local enlargement is used... This achieves adaptive reduction of the threshold spacing, thereby forming a local encryption layer.

[0076] After obtaining the threshold, isosurface families are extracted along the trans-wall harmonic function. The voxel mesh uses the isosurface cube method, and the tetrahedral mesh uses the triangulation method, outputting closed patches with consistent normals. To avoid jagged edges and hanging triangles, in-plane smoothing and short-side collapse are performed first, followed by shape-preserving constraints at high curvature locations to maintain geometric details. The actual area of ​​each layer's patches is calculated and compared with the target area. If the deviation exceeds the limit, the process returns to the threshold calculation step to fine-tune adjacent thresholds until the layer area deviation meets the preset requirements.

[0077] To generate polygonal surface mesh units aligned with the anatomical direction, local coordinates are defined using the gradient directions of three sets of harmonic functions. On the isosurface of each layer, isolines corresponding to the direction from the interventricular septum to the lateral wall are selected as meridians, and isolines corresponding to the direction from the apex to the base are selected as parallels. The intersection of meridians and parallels forms corner points, and connecting adjacent corner points yields quadrilateral or oligogonal surface units with better regularity. For locally geometrically complex regions, such as near the valve annulus and outflow tract, a denser spacing of meridians and parallels is used, and the distortion and angle range of the surface units are limited; if abnormally elongated surface units appear, local re-partitioning is triggered.

[0078] The implementation details are as follows: The curvature component of the weight field is derived from the second-order difference of the transwall harmonic function; the region indicator function is derived from the labeled lobe rings and outflow duct masks, and undergoes morphological dilation to ensure coverage edges; the weight coefficients are automatically calibrated using a small number of mesh quality samples. The cumulative volume function is evaluated more quickly using prefix sum voxel integration; the threshold solution uses multi-point parallel bisection to ensure global monotonicity. The triangulation stage of the isosurface records normal consistency and performs point-to-surface nearest-point registration between layers to improve the stability of inter-layer correspondence. The generation of latitude and longitude lines combines in-surface level set tracking with terrain line tracking, automatically reducing the step size at abrupt normal changes to avoid exceeding boundaries.

[0079] Quality control includes three types of checks. First, area deviation check: the relative error between the area of ​​each layer and the target area must be within a threshold; otherwise, it returns to the threshold fine-tuning stage for iterative correction. Second, direction consistency check: the angles between the tangents of the meridians and parallels and the gradient directions of the two sets of harmonic functions should be less than the set upper limit; if the limit is exceeded, the line density and line direction are adjusted locally. Third, mesh shape check: the minimum angle, maximum angle, side length ratio, and normal fluctuation of the surface elements must all be qualified simultaneously; unqualified surface elements trigger local re-division or shape-preserving smoothing. The final output is a polygonal surface mesh with controlled layer area, local densification in key areas, and meridians and parallels consistent with the anatomical direction, which can directly enter the subsequent closed-loop optimization and application solution.

[0080] Using homeomorphic deformation and three sets of harmonic function fields as joint unknowns, a target including pixel consistency, physical consistency and topological consistency is constructed. The polygonal surface mesh is updated and reconstructed synchronously. The density of the inducer is adjusted based on the layer area deviation feedback until convergence, and the three-dimensional polygonal surface mesh of the heart is output.

[0081] This embodiment uses homeomorphic deformation and three sets of harmonic fields as joint unknowns (joint unknowns refer to the set of variables solved simultaneously in the same optimization process, including homeomorphic deformation mapping and three sets of harmonic fields), and adopts a two-layer closed-loop optimization. In the inner layer, under the condition of a fixed layer threshold and the current mesh, the homeomorphic deformation and three sets of harmonic fields are updated synchronously: first, pixel residuals are generated using differentiable reprojection and backpropagated to voxel coordinates; then, physical residual gradients are generated from the physical adjoint solution, and the two are weighted and mapped to the velocity field and harmonic field; in the inaccessible zone, the normal component is zeroed out, and only the tangential component is updated; local oscillations are suppressed using step size control and anisotropic smoothing; after each update, anchor point position correction and Jacobian determinant lower bound checks are performed, and if not satisfied, backtracking and reducing the step size are performed. Subsequently, isosurfaces are extracted according to the latest three sets of harmonic fields, and polygonal surface meshes are reconstructed. Local repartitioning and shape-preserving smoothing are performed on small-angle and slender surface elements to maintain manifold and anatomical alignment.

[0082] The outer layer, used for area constraint feedback and induced volume density reweighting, is calculated as follows: ; , in, For the first Relative deviation of layer area For the first Measured area of ​​each floor For the first Target area of ​​the layer For the density of the inducible body, The reweighted inducible body density. This is the step size coefficient. After weighting, the mass-conserving transport is re-executed to update the transport map, providing a consistent data distribution basis for the next round of inner-layer synchronous updates.

[0083] In practice, multi-scale stepwise optimization is used to expand the convergence domain; during the isosurface reconstruction stage, the normal direction is unified and the correspondence between points and surfaces between layers is maintained to reduce interlayer misalignment; directional smoothing is used in the vicinity of the valve annulus and outflow tract to limit normal stretching; quality control is based on the target descent rate, maximum anchor point deviation, minimum Jacobian determinant, relative deviation of layer area, and mesh quality indicators, all of which must be met simultaneously for convergence to occur. The final output is the converged homeomorphic deformation, three sets of harmonic function fields, and the three-dimensional polygonal surface mesh of the heart, along with its quality control report.

[0084] Synchronous updates involve constructing a block linear system with Karush-Kuhn-Tucker constraints using homeomorphic deformation and three sets of harmonic function fields as joint unknowns. The constraints include pixel consistency constraints, physical consistency constraints, and topological consistency constraints. An iterative linear solver is used to solve the block linear system to achieve synchronous updates and reconstruction of the polygonal surface mesh.

[0085] This embodiment illustrates the synchronous update and reconstruction process. Homeomorphic deformation and three sets of harmonic function fields are used as joint unknowns to construct a KKT (Karush-Kuhn-Tucker) block linear system containing pixel consistency, physical consistency, and topological consistency. An iterative linear solver is used to achieve a closed loop of solving and reconstructing in one step. The calculation expressions are as follows: The Hessian matrix is ​​the Gaussian-Newton approximation for homeomorphic deformation. For the Gaussian-Newton approximation Hessian matrix of three sets of harmonic function fields, The Jacobian matrix constrained by pixel consistency. The Jacobian matrix is ​​a physical consistency constraint. The Jacobian matrix is ​​a topological consistency constraint. For homeomorphic deformation increment, For the three sets of harmonic function field increments, For pixel consistency Lagrange multipliers, For physically consistent Lagrange multipliers, For topologically consistent Lagrange multipliers, The first gradient of homeomorphic deformation. The first gradients of the three harmonic function fields are: For reprojection residuals, For physical residuals, This represents the topological residual. The construction of the above pixel and physical terms is described in the aforementioned joint formula and will not be repeated here.

[0086] The implementation process is as follows: Assemble the blocks of matrices and vectors at the current iteration point based on automatic differentiation; solve using MINRES or GMRES for symmetric indeterminate problems, and use a block diagonal-Shure complement preconditioner to reduce the number of iterations; obtain the solution. and Then, the homeomorphic deformation and harmonic function fields are updated by line search; the normal component is set to zero in the uncrossable zone, and only the tangential update is retained; a threshold barrier check is performed on the Jacobian determinant, triggering backtracking and step size reduction until the positive value constraint is met; local least squares orthogonalization is performed on the three sets of harmonic function fields and the boundary values ​​are written back to maintain consistency with the dissection direction.

[0087] Immediately after each linear solution, the polygonal surface mesh is reconstructed: isosurfaces are extracted based on the latest harmonic field, corner points are constructed along the isosurfaces in three directions, and surface elements are generated; edge collapse is performed on short sides, and shape-preserving smoothing is applied to acute-angle surface elements; only tangential smoothing is allowed in the neighborhood of the lobe ring and outflow channel to protect the anatomical curves. Subsequently, the layer area is calculated and compared with the target area, and the area deviation curve is recorded as the input for outer layer feedback; if the mesh quality or area deviation does not meet the standard, the next round of assembly-solution-update-reconstruction continues while maintaining the validity of constraints. This process synchronously couples geometric and field variables within a single linear system, reducing the risk of mismatch in step-by-step optimization and ensuring that pixel, physical, and topological consistency converge together under a unified framework.

[0088] Adjusting the induced volume density based on layer area deviation feedback involves converting the layer area deviation into a density reweighting factor, reweighting the induced volume density, and using it as the target density in the next mass conservation transmission. This process is repeated synchronously until the polygonal surface mesh meets the preset mesh quality and topology consistency requirements.

[0089] This embodiment addresses closed-loop correction of layer area deviation by providing an engineering workflow that transforms the layer area deviation into inducible volume density, reweights it, and drives the next round of mass-conserving transport. The inputs are the polygonal surface mesh generated in the previous round, the layer threshold, the inducible volume density, and the patient-specific cardiac 3D domain; the outputs are the reweighted inducible volume density, the updated transport map, and the reconstructed polygonal surface mesh, until the mesh quality and topology consistency requirements are met.

[0090] First, the surface area is statistically analyzed within each transwall layer. Normal unification and orifice repair are performed on the isosurfaces, and the area is calculated using surface triangulation results. Independent statistics are provided for the vicinity of the lobe ring and outflow channel to avoid interference from local geometry on the overall evaluation. After obtaining the measured area layer by layer, the relative deviation of the layer area is calculated according to the aforementioned "layer area feedback expression" (described above), and this deviation is used to generate the density reweighting factor for the layer. To suppress numerical oscillations, truncation and one-dimensional monotonic filtering are applied to the reweighting factor to ensure continuous factor changes between adjacent layers; narrow-band smoothing is used inside the layer boundaries to avoid abrupt changes in induced volume density; and mass normalization is performed on the global induced volume density to ensure that the total mass remains constant and the voxel values ​​are non-negative. The meanings of relevant symbols and the definition of the step size coefficient γ are explained in the descriptions below the aforementioned expressions and will not be repeated here.

[0091] After reweighting, the mass conservation transport update step is invoked, replacing the target density with the reweighted induced volume density. The solution process is consistent with the "transport solution expression" (described earlier), using the previous round's proportional vector as initialization and maintaining the sparse neighborhood structure of the core at the layer boundary to stabilize the transport plan. After obtaining the new transport map, it is used as the initial value for homeomorphic deformation and the three sets of harmonic function fields to drive a synchronous update. The update method is consistent with the aforementioned KKT block linear system solution, and the algorithm details will not be repeated. After the update, isosurfaces are extracted according to the latest harmonic function fields, the polygonal surface mesh is reconstructed, and three types of quality control indicators are calculated: relative deviation of layer area, mesh shape index, and topology consistency index. The mesh shape index includes minimum angle, maximum angle, side length ratio, and surface normal fluctuation; the topology consistency index includes the lower bound of the Jacobian determinant, span count, and manifold check.

[0092] If any indicator fails to meet the standard, the next outer layer cycle begins: the step size coefficient γ and the upper limit of the reweighting factor are automatically adjusted based on the new relative deviation of the layer area, prioritizing the increase of correction intensity for layers with stable deviation signs and large amplitudes; stronger narrowband smoothing and smaller step sizes are used for high curvature regions and the vicinity of the valve annulus and outflow tract to avoid local distortion caused by over-correction. The cycle terminates when the relative deviation of the entire layer area does not exceed the threshold, all mesh shape indicators are qualified, and the topology consistency check is passed. The final output includes a three-dimensional polygonal surface mesh of the heart, the corresponding layer threshold, the reweighted inducible body density and transport mapping, and a quality control report covering the above indicators, which can be directly used for subsequent numerical analysis and visualization.

[0093] like Figure 2 As shown, a cardiac 3D polygon mesh construction system is used to implement the aforementioned cardiac 3D polygon mesh construction method. The system includes: The image modeling module is used to acquire and register medical images, construct inducible body density and signed distance fields based on boundary features, and fit the anatomical anchor point set. It is recommended that the image modeling module utilize a workstation or edge server with a multi-core general-purpose processor and a single / dual-processor graphics processor. The processor handles DICOM reception, decoding, and I / O scheduling, while the graphics processor performs parallel tasks such as voxel-level filtering, gradient and curvature calculations. A large capacity of memory is configured for resident 3D volume data and multi-scale pyramid caching, with NVMe SSDs used for high-speed intermediate results and slice caching, and RAID 1 / 10 configuration to ensure data reliability. The network side is equipped with 10GbE or higher Ethernet ports, supporting DICOM over TLS and direct connection to PACS; optional hardware time synchronization (e.g., IEEE 1588) is available for multi-device collaboration. A medical-calibrated monitor is used for display. The entire system must have hot-swappable drive bays and redundant power supplies to ensure stability for continuous image library writing and long-term operation.

[0094] The transport deformation module is used to calculate mass-conserving transport between the template volume density and the inducer volume density to obtain the transport map, calculate homeomorphic deformation and constrain it at the anatomical anchor set, and combine pixel consistency and physical consistency to determine the patient-specific cardiac 3D domain. Numerical optimization of mass-conserving transport, differential reprojection, and homeomorphic deformation is sensitive to video memory and bandwidth, so a high-bandwidth graphics processor (HBM architecture) with multi-GPU support is recommended. This module is used for kernel matrix operations, scaling iterations, and velocity field exponential mapping. The host is equipped with a high-frequency multi-core processor responsible for sparse adjacency construction, boundary condition orchestration, and constraint checking. ≥256GB of memory is recommended to accommodate sparse matrix and multi-resolution volume data. Multiple NVMe drives are used in parallel on the system disk and Scratch disk to improve random read / write speeds. The machine provides a PCIe Gen4 / Gen5 bus and ample power supply and cooling. A 25GbE uplink network is provided to support data backhaul between remote PACS and compute nodes.

[0095] The metric meshing module is used to construct pullback metrics within the patient-specific 3D cardiac domain, solve for three sets of harmonic function fields aligned with the anatomical direction, extract isosurfaces based on co-area constraints, and generate polygonal surface meshes through their intersection. Pullback metric assembly, finite element discretization, and harmonic function solving all carry significant computational and memory access pressures; therefore, a large-memory, multi-core server is recommended as the primary platform. The CPU's vectorization capabilities should be utilized to complete the sparse linear system assembly and algebraic multimesh preconditioning. A graphics processor is used for parallel gradient field and isosurface triangulation, surface tracing, and other high-concurrency tasks. To improve mesh quality inspection and geometric operation efficiency, a high-speed NVMe repository is deployed as a local repository for mesh versions and checkpoints, and a UPS and ECC memory are configured to ensure power and data reliability during long-term iterations. If necessary, 3D input devices and high-resolution display terminals are provided for manual review and interactive correction.

[0096] The closed-loop optimization module constructs an objective network encompassing pixel consistency, physical consistency, and topological consistency using homeomorphic deformation and three sets of harmonic function fields as joint unknowns. It synchronously updates and reconstructs the polygonal surface network, adjusting the inducible body density based on layer area deviation feedback until convergence, outputting a three-dimensional polygonal surface network of the heart. Solving the KKT block linear system with joint unknowns and multiple rounds of outer-layer feedback requires stable computational orchestration and high-throughput I / O. It is recommended to use a computing node (CPU + multiple GPUs) with job scheduling and fault tolerance as the optimization host. GPUs handle Jacobi / gradient batch calculations and synchronous updates, while CPUs handle KKT assembly, line search, and quality control indicator statistics. On the storage side, high-speed local NVMe combined with centralized NAS / SAN is used for checkpoint and quality report archiving; the network is at least 25GbE to ensure low-latency data exchange between multiple modules. To ensure availability in the medical environment, the entire system is equipped with redundant power supplies, rack-level cooling and noise control, and reserved hardware security modules and encryption acceleration to support full lifecycle data encryption and access auditing, ultimately achieving a hardware closed-loop execution of multiple rounds of update-reconstruction-feedback.

[0097] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects.

[0098] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A method for constructing a three-dimensional polygonal mesh of the heart, characterized in that, Includes the following steps: Acquire and register medical images, construct inducible body density and signed distance field based on boundary features, and fit the set of anatomical anchor points; Mass conservation transport is calculated between template body density and inducer body density to obtain transport mapping, homeomorphic deformation is calculated and constrained at the anatomical anchor set, and the patient-specific cardiac three-dimensional domain is determined by combining pixel consistency and physical consistency. A pullback metric is constructed in the three-dimensional domain of the patient's specific heart. Three sets of harmonic function fields are solved and aligned with the anatomical direction. Isosurfaces are extracted based on co-area constraints and their intersections are used to generate a polygonal surface mesh. Using homeomorphic deformation and three sets of harmonic function fields as joint unknowns, a target including pixel consistency, physical consistency and topological consistency is constructed. The polygonal surface mesh is updated and reconstructed synchronously. The density of the inducer is adjusted based on the layer area deviation feedback until convergence, and the three-dimensional polygonal surface mesh of the heart is output.

2. The method according to claim 1, characterized in that, Acquiring and registering medical images includes spatial registration and intensity normalization, extracting boundary features based on gradient structure tensor and curvature estimation, constructing inducible body density from boundary features through kernel smoothing, calculating the signed distance field from the topologically purified initial boundary using the fast travel method, and the set of anatomical anchor points includes the valve annulus closure curve, the apex of the heart, the opening loop of the great vessels, and the master curve of the interventricular septum.

3. The method according to claim 1, characterized in that, The calculation of mass-conserving transport between the template volume density and the induced volume density involves solving for the optimal transport using entropy regularization, with the cost function being the squared Euclidean distance.

4. The method according to claim 3, characterized in that, The calculation of homeomorphic deformation includes calculating homeomorphic deformation based on the exponential mapping of the static velocity field, applying positional constraints at the set of anatomical anchor points, setting non-crossable boundary constraints in the atrioventricular septal region, and setting positive constraints on the Jacobian determinant to maintain the mapping as homeomorphic.

5. The method according to claim 1, characterized in that, Pixel consistency is calculated from the reprojection error of differentiable rendering, while physical consistency is calculated from the residuals of the forward and adjoint solutions of the wave equation. The gradients of the two are backpropagated along the mapping relationship of homeomorphism to update the homeomorphism and determine the patient-specific cardiac three-dimensional domain.

6. The method according to claim 1, characterized in that, The construction of the pull-back metric involves pulling back the Euclidean metric through homeomorphic deformation to obtain the metric tensor. Three sets of harmonic function fields are subjected to Dirichlet boundary conditions at the apex and base, Dirichlet boundary conditions at the interventricular septum and Neumann boundary conditions at the lateral wall, Dirichlet boundary conditions at the endocardium and epicardium, and Neumann boundary conditions at the opening of the great vessels.

7. The method according to claim 6, characterized in that, The extraction of isosurfaces based on co-area constraints includes constructing weights based on curvature and anatomical region indicator functions. According to the weight allocation level, the weights are increased in the valve ring and outflow channel regions to form local densification, and polygonal surface mesh units aligned with the anatomical direction are generated on the intersection line of the isosurfaces.

8. The method according to claim 1, characterized in that, Synchronous updates involve constructing a block linear system with Karush-Kuhn-Tucker constraints using homeomorphic deformation and three sets of harmonic function fields as joint unknowns. The constraints include pixel consistency constraints, physical consistency constraints, and topological consistency constraints. An iterative linear solver is used to solve the block linear system to achieve synchronous updates and reconstruction of the polygonal surface mesh.

9. The method according to claim 8, characterized in that, Adjusting the induced volume density based on layer area deviation feedback involves converting the layer area deviation into a density reweighting factor, reweighting the induced volume density, and using it as the target density in the next mass conservation transmission. This process is repeated synchronously until the polygonal surface mesh meets the preset mesh quality and topology consistency requirements.

10. A cardiac three-dimensional polygonal mesh construction system, used to implement the cardiac three-dimensional polygonal mesh construction method according to any one of claims 1 to 9, characterized in that, The system includes: The image modeling module is used to acquire and register medical images, construct inducible body density and signed distance fields based on boundary features, and fit the set of anatomical anchor points. The transport deformation module is used to calculate the mass conservation transport between the template body density and the inducer body density to obtain the transport mapping, calculate the homeomorphic deformation and constrain it at the anatomical anchor point set, and combine pixel consistency and physical consistency to determine the patient-specific cardiac three-dimensional domain. The metric meshing module is used to construct pullback metrics in the three-dimensional domain of the patient's specific heart, solve three sets of harmonic function fields and align them with the anatomical direction, extract isosurfaces based on co-area constraints and generate polygonal surface meshes by their intersection; The closed-loop optimization module is used to construct a target that includes pixel consistency, physical consistency and topological consistency using homeomorphic deformation and three sets of harmonic function fields as joint unknowns. It synchronously updates and reconstructs the polygonal surface mesh, adjusts the inducing body density based on layer area deviation feedback until convergence, and outputs a three-dimensional polygonal surface mesh of the heart.

Images