Heart failure follow-up data visualization processing method based on multi-modal image fusion

Abstract

The present application relates to the technical field of image fusion, and discloses a heart failure follow-up data visualization processing method based on multi-modal image fusion, comprising the following steps: constructing a left ventricular fiber direction field and extracting fusion scalar data along the fiber path, establishing an isochronous scalar complex quantification model for describing the change pattern of activation time in the fiber direction, and taking the model as a continuous control variable to participate in subsequent filtering, reconstruction and visualization coding processes. The method completes the structured arrangement of the scalar data along the integral path of the fiber, so that the fusion data is more consistent with the propagation law of the myocardial tissue when outputting the polar coordinate graph and the three-dimensional ventricular coloring graph, and can visually inhibit and prompt the local non-continuous fluctuation.

Description

Technical Field

[0001] This invention relates to the field of image fusion technology, and more specifically, to a method for visualizing heart failure follow-up data based on multimodal image fusion. Background Technology

[0002] In the diagnosis, prognostic assessment, and long-term follow-up of heart failure, obtaining information on myocardial deformation and mechanical activation using cardiac magnetic resonance imaging (CMRcine, tagging, DENSE), ultrasound speckle tracking, and other medical imaging techniques has become routine. Different imaging modalities vary in spatial resolution, temporal resolution, signal-to-noise ratio, and sensitivity. Therefore, in clinical practice, it is often necessary to integrate multiple imaging results into a more intuitive and unified analysis to support the observation of key physiological changes such as myocardial contractile dyssynchrony, delayed mechanical activation regions, and myocardial remodeling progression. Currently widely used visualization methods mainly include: AHA-based segmented polar coordinate plots (bull's-eye plots), three-dimensional ventricular surface shading maps, and volumetric two-dimensional unfolded plots. These can display strain distribution, activation time delay distribution, and their follow-up changes, enabling physicians to quickly assess cardiac function and track the evolution of myocardial dynamic patterns over time.

[0003] While the methods described above can integrate image information from different sources, existing visualization workflows generally only focus on the distribution of the fused scalar data in two-dimensional or three-dimensional space, without fully considering the inherent structural propagation patterns within the myocardium. For example, myocardial tissue has a spiral arrangement of muscle fibers from the endocardium to the adventitia, and mechanoactivation actions typically proceed continuously along these fiber directions or fiber-dominated paths. However, when data from different modalities are affected by noise perturbations or response differences during spatial reconstruction, temporal registration, and interpolation fusion, the scalar data may exhibit irregular fluctuations along the myocardial fiber directions, or even form repetitive local oscillations. Because existing systems lack a processing workflow to analyze along the actual conduction paths of the myocardium, these perturbations are visually averaged in the final bull's-eye map or 3D shaded map, appearing smooth and thus masking potential variations within the data.

[0004] Specifically, the propagation of mechanical forces within the myocardium does not occur arbitrarily in Euclidean space, but is strictly constrained by natural histological laws such as fiber orientation, interlaminar rotation structure, and ventricular wall thickness. When the activation time or strain distribution after fusion violates these histological laws and produces local distortions, if the visualization process simply colors the scalar data based on geometric mapping, it cannot express the degree of abnormality of these local changes in a mechanical sense. Summary of the Invention

[0005] This invention provides a method for visualizing heart failure follow-up data based on multimodal image fusion, which solves the technical problems mentioned in the background art.

[0006] This invention provides a method for visualizing heart failure follow-up data based on multimodal image fusion, including:

[0007] A geometric model of the left ventricle and a three-dimensional myocardial fiber architecture model were constructed based on multimodal follow-up images, and initial full left ventricular myocardial biomechanical activation distribution data were generated by fusing them.

[0008] A virtual conduction path is constructed along the three-dimensional myocardial fiber architecture model, and the distribution of fiber conduction topological disorder, which characterizes the degree of non-physical re-entry of the activation wave along the fiber, is calculated.

[0009] Based on the distribution of fiber conduction topology disorder, the filtering parameters are adjusted, and adaptive smoothing is performed only along the virtual conduction path on the whole left ventricular myocardial biomechanical activation distribution data to generate myocardial activation distribution data corrected by fiber structure constraints.

[0010] The corrected residual conduction disorder distribution is calculated and mapped in a multidimensional visual manner with the myocardial activation distribution data corrected by fiber structure constraints, generating a heart failure follow-up visualization image that simultaneously presents the activation sequence and structural reliability on different visual channels.

[0011] The beneficial effects of this invention include: by introducing isochronous scalar complexity analysis along the myocardial fiber path, the activation time or strain information after multimodal fusion is continuously corrected and structurally visualized according to the actual tissue structure of the myocardium. This allows the visualization results to no longer rely solely on color mapping itself, but to reflect the variation patterns and reliability of scalar data along the fiber direction, thereby significantly improving the identifiability and controllability of local abnormal fluctuations during follow-up. The final output polar coordinate map and three-dimensional ventricular shading map are visually closer to the myocardial dynamics propagation mechanism, effectively avoiding misinterpretation of changes caused by multimodal fusion errors. This provides a more robust, intuitive, and structurally consistent presentation method for identifying myocardial dynamics trends and tracking treatment responses in long-term clinical follow-up. Attached Figure Description

[0012] Figure 1 This is a flowchart of the heart failure follow-up data visualization processing method based on multimodal image fusion according to the present invention;

[0013] Figure 2 This is a schematic diagram of the fiber field and integral trajectory of the present invention;

[0014] Figure 3 This is a schematic diagram illustrating the principle of calculating the isochronous scalar Morse complexity along the myocardial fiber path according to the present invention. Detailed Implementation

[0015] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0016] like Figure 1 As shown, the method for visualizing heart failure follow-up data based on multimodal image fusion includes:

[0017] A geometric model of the left ventricle and a three-dimensional myocardial fiber architecture model were constructed based on multimodal follow-up images, and initial full left ventricular myocardial biomechanical activation distribution data were generated by fusing them.

[0018] A virtual conduction path is constructed along the three-dimensional myocardial fiber architecture model, and the distribution of fiber conduction topological disorder, which characterizes the degree of non-physical re-entry of the activation wave along the fiber, is calculated.

[0019] Based on the distribution of fiber conduction topology disorder, the filtering parameters are adjusted, and adaptive smoothing is performed only along the virtual conduction path on the whole left ventricular myocardial biomechanical activation distribution data to generate myocardial activation distribution data corrected by fiber structure constraints.

[0020] The corrected residual conduction disorder distribution is calculated and mapped in a multidimensional visual manner with the myocardial activation distribution data corrected by fiber structure constraints, generating a heart failure follow-up visualization image that simultaneously presents the activation sequence and structural reliability on different visual channels.

[0021] In a preferred embodiment, a left ventricular geometric model and a three-dimensional myocardial fiber architecture model are constructed based on multimodal follow-up images, including:

[0022] Select the left ventricular reference geometry grid in multimodal follow-up images. Calculate the value of each point within the grid. Normalized trans-wall depth The endometrium outer membrane The normalized transmembrane depth is defined as the linearly normalized distance from the inner membrane to the outer membrane;

[0023] Calculate the local helix angle at each point using the following formula. :

[0024] ;

[0025] in, The intima helix angle, The outer membrane helix angle;

[0026] The three-dimensional myocardial fiber architecture model was synthesized and normalized according to the following formula. :

[0027] ;

[0028] in, Let be the unit vector along the circumference in the local cylindrical coordinate system. Let be the unit vector along the vertical axis of the local cylindrical coordinate system.

[0029] It should be noted that existing technologies typically process multimodal data directly in Cartesian coordinates or simple polar coordinates, neglecting the fact that myocardial electromechanical waves propagate along the myocardial fibers. Therefore, by constructing a three-dimensional myocardial fiber architecture model, we can ensure that the isochronous Morse complexity calculated along the fiber path in subsequent calculations accurately reflects the topological structure of myocardial mechanical conduction, rather than simply spatial geometric features.

[0030] Preferably, a mid-diastolic image from a cardiac magnetic resonance (CMR) cinema sequence is selected as the reference image, as the left ventricle appears relatively full and motion artifacts are minimal at this time. Using a U-Net-based deep learning segmentation network or a conventional level set segmentation algorithm, the endocardial and epicardial boundaries of the left ventricle are extracted, and a left ventricular reference geometric mesh composed of tetrahedral or hexahedral units is constructed. For any point in space within the grid First, calculate its normalized cross-wall depth. Specifically, this depth is defined as the shortest Euclidean distance from this point to the inner membrane surface. The local myocardial wall thickness at that point The ratio, that is ,in It can be approximated as the sum of the shortest distance from the point to the inner membrane and the shortest distance to the outer membrane. It is a dimensionless scalar whose range of values ​​is restricted to a closed interval. Inside, among which Indicates the inner membrane surface. The term represents the surface of the adventitia, and this parameter is used to parameterize the thickness direction of the myocardial wall. After obtaining the normalized transmural depth, a modified Streeter helical fiber model is used to analytically construct the three-dimensional myocardial fiber orientation. Based on prior anatomical knowledge, the helix angle of myocardial fibers exhibits a continuous linear trend from a positive angle in the endocardium to a negative angle in the adventitia. In specific calculations, the endocardial helix angle parameter is preset. for outer membrane helix angle parameter for These two parameters are empirical constants derived from a large number of in vitro cardiac tissue measurements. For any point within the grid... Using linear interpolation formula Calculate the local helix angle at this point. . To normalize the cross-wall depth, and The preset angular boundary values ​​are then used. Subsequently, the scalar helix angle is transformed into a unit vector field in three-dimensional space. In the local cylindrical coordinate system of the left ventricle, a local unit vector along the longitudinal axis is defined. (Pointing to the base of the left ventricle) and the unit vector in the local circumferential direction (In a counter-clockwise direction along the tangential axis). Three-dimensional myocardial fiber architecture model. Through formula Perform synthesis. In this formula, The term represents the projection component of the fiber vector in the circumferential direction. The term represents the projection component of the fiber vector along the vertical axis. Since The generated vector Satisfy the unit modulus constraint, i.e. This allows for the construction of a continuous, smooth, and differentiable vector field on a discrete geometric grid. This vector field not only reflects the macroscopic geometry of the left ventricle but also embeds the microscopic spiral arrangement of myocardial fibers.

[0031] In some possible embodiments, the parameters and models in the above construction process can be adaptively adjusted to cover more technical variations and adapt to the pathological characteristics of different patients. For example, in patients with dilated cardiomyopathy (DCM), the normalized transmural depth is adjusted to account for the potential for ventricular remodeling to alter fiber arrangement. The calculation can be performed using the solution of the Laplace equation instead of a simple distance ratio. That is, the potential energy at the inner membrane boundary is set to 0, and the potential energy at the outer membrane boundary is set to 1. Solving for the solution... The resulting potential field is used as the depth field, making this method smoother and more robust when dealing with regions of significant myocardial wall curvature variation. Furthermore, the parameters in the helix angle model... and It can also be corrected based on the patient's personalized diffusion tensor magnetic resonance imaging (DT-MRI) data, or set to... to Full coverage to accommodate different myocardial structures. Simultaneously, the fiber vector synthesis formula can also incorporate a transverse angle component, making... In the radial direction It also has a component, thus enabling the construction of a more complex generalized myocardial fiber model.

[0032] In a preferred embodiment, initial full left ventricular myocardial biomechanical activation distribution data are generated by fusing the left ventricular geometric model and the three-dimensional myocardial fiber architecture model, including:

[0033] Will Data from each imaging modality is registered to the left ventricular geometric model, and mechanical activation time data for each modality are extracted. ,in Indicates the first One mode;

[0034] The distribution data of total left ventricular myocardial biomechanical activation were calculated using the following formula. :

[0035] ;

[0036] in, For the first Each modality has a preset fixed fusion weight, and satisfies as well as .

[0037] It should be noted that existing technologies, when processing multimodal data (such as cardiac MRI and echocardiography), often introduce complex consistency determination or optimal selection logic, i.e., when different modal data conflict, a choice is made between two or the former is eliminated. This nonlinear logic processing can lead to numerical jumps or holes in the generated activation map in space, interfering with subsequent differential calculations and topological analysis along the fiber direction. Therefore, by adopting a fixed-weight linear fusion strategy, the aim is to accept input information from all modalities without discrimination.

[0038] Preferably, the fusion process first involves the spatial and temporal registration of the multimodal images. For example, the input image modalities include cardiac magnetic resonance displacement-coded imaging (CMR-DENSE, denoted as modality). ) and three-dimensional speckle tracking echocardiography (3D-STE, denoted as modal) First, using the spatial positioning information in the DICOM header file and a rigid registration algorithm based on mutual information, the coordinate systems of each mode are aligned to the left ventricular geometric model. For any node on the geometric model mesh... The trilinear interpolation method is used to extract the corresponding mechanical activation time values ​​from the registered original data volumes of each mode, thereby forming a single-mode activation distribution field. and Here, The unit is milliseconds (ms), representing the time delay from the onset of the QRS wave on an electrocardiogram to the point at which the maximum principal strain (or maximum systolic rate) is reached. If there are differences in heart rate between different modalities, the time values ​​need to be normalized to the same cardiac cycle scale beforehand based on the RR interval. Subsequently, a linear weighted summation calculation is performed. The calculation formula is as follows: In this formula, This refers to the initial full left ventricular myocardial biomechanical activation distribution data to be generated, which represents the estimated local activation time after fusion. and These are preset fixed fusion weights. Specifically, the parameters... and The value of is determined based on prior knowledge of the spatial resolution and signal-to-noise ratio of each modality. For example, given that CMR-DENSE is generally more accurate than ultrasound in three-dimensional myocardial displacement measurement and is not limited by the acoustic window, is set to . ,set up Both parameters are non-negative constants and strictly satisfy the normalization constraint. .

[0039] In some possible embodiments, cardiac magnetic resonance labeling (CMR-Tagging, denoted as...) is included. In a three-modal fusion scenario, the formula is expanded to: In this case, the weights can be allocated based on the dimensional characteristics of the modality. For example, since both DENSE and Tagging are three-dimensional acquisitions, their weights can be set higher (e.g., 0.4 each), while ultrasound, being an auxiliary modality, has a lower weight (e.g., 0.2). In another variation, considering the differences in imaging quality between different modalities in different ventricular segments (e.g., ultrasound imaging is better at the apex, while MRI is better at the base), the fixed weights can be implemented as a spatially location-dependent preset weight field. That is, for a specific anatomical region (such as the apical segment of the AHA), W is pre-set. echo Greater than W cmr The opposite is true for the basal segment.

[0040] In a preferred embodiment, constructing a virtual conduction path along the three-dimensional myocardial fiber architecture model includes:

[0041] Uniformly selected within the left ventricular myocardial wall of the left ventricular geometric model Representative seed points ( (index of the seed point).

[0042] For each representative seed point Along the three-dimensional myocardial fiber architecture model The virtual conduction path is obtained by bidirectional extension using a numerical integration method with a fixed step size. It satisfies the following differential equation:

[0043] ;

[0044] in, Let be the arc length parameter along the path, and the fixed step size for numerical integration is . Stop extending beyond the boundaries of the myocardium or when the maximum number of steps is reached;

[0045] The numerical values ​​of the whole left ventricular myocardial biomechanical activation distribution data are sampled at fixed intervals along each of the virtual conduction paths to obtain a one-dimensional activation function distributed along the path arc length. :

[0046] ;

[0047] in, Data on the distribution of mechanical activation in the entire left ventricular myocardium. Indicates the first The arc length coordinates of each sampling point.

[0048] It should be noted that existing techniques typically calculate the magnitude or direction of the gradient directly in three-dimensional Euclidean space. However, this ignores the unique characteristics of the myocardial medium. Electroexcitation and mechanical contraction waves do not diffuse spherically in an isotropic medium, but are strongly constrained by helical myocardial fibers, primarily propagating along the fiber direction and its orthogonal lamellar directions. Therefore, constructing a virtual conduction path through numerical integration essentially simulates the potential propagation trajectory of the physical wavefront. Analyzing the changes in activation time along this trajectory can shield against spatial heterogeneity interference perpendicular to the fiber direction, thus transforming the three-dimensional topological problem into a more easily analyzed one-dimensional signal monotonicity problem.

[0049] Preferably, the process of constructing the virtual conduction pathway begins with the selection of seed points. To ensure that the pathway can uniformly cover all anatomical regions of the left ventricular wall, voxel sampling is performed based on the aforementioned constructed left ventricular geometric grid. Specifically, sampling density parameters are set. (e.g., 2.0 mm) A uniform three-dimensional lattice is generated within the region of interest (ROI) of the myocardium surrounded by the endocardium and adventitia, and each point in this lattice is used as a fiber origin. ,in , This represents the total number of seed points. Then, for each starting point... In the three-dimensional myocardial fiber architecture model Two-way numerical integration is performed. The integration process follows the first-order ordinary differential equation: The initial conditions are To ensure the geometric accuracy of the integral trajectory and suppress accumulated errors, the fourth-order Runge-Kutta method (RK4) is preferred for iterative solution. The single-step iteration formula is as follows: ,in The preset integration step size is typically set to the minimum size of the image voxels. to Times (e.g., 1.0 mm); to Based on the vector field at different intermediate points within the current step size The slope of the calculated tangent vector. Integration is performed simultaneously along the positive (+f) and negative (-f) directions until the trajectory touches the boundaries of the endocardium, epicardium, basal plane, or apex of the left ventricular geometry, or the trajectory length exceeds the preset maximum arc length L. max (For example, 150 mm), at which point the integration ends, yielding a complete space curve. Finally, data sampling is performed along the path. For each discrete integration point on the curve... The trilinear interpolation algorithm was used to analyze the mechanical activation distribution data of the entire left ventricular myocardium. The corresponding activation time values ​​are extracted from the voxel mesh, thereby transforming the three-dimensional scalar field. Transform into a parameterized one-dimensional activation function .

[0050] In some possible implementations, the Euler method, which is more computationally efficient but slightly less accurate, can be used instead of the Runge-Kutta method, especially for rapid initial screening of large-scale clinical data. Specifically, the iterative formula simplifies to... Although its local truncation error is relatively large, the fiber orientation field is relatively smooth and the integration step size is large. Even with sufficiently small seed point settings (e.g., 0.2 mm), paths that meet the requirements of topology analysis can still be generated. Furthermore, regarding seed point selection, in addition to uniform grid sampling, a center-point sampling strategy based on the AHA17 segment model can be used. This involves placing seed points at the centroid of each myocardial segment and within a small surrounding area. This generates fewer paths but with clear anatomical representativeness, significantly reducing computational load. In the sampling stage, if the input multimodal data has low resolution, cubic spline interpolation can be used instead of linear interpolation to obtain a one-dimensional function. Continuity of higher-order derivatives.

[0051] In a preferred embodiment, calculating the fiber-borne topological disorder distribution characterizing the degree of non-physical back-reflection of the active wave propagating along the fiber includes:

[0052] Calculating the one-dimensional activation function using the finite difference method The derivative with respect to arc length And determine the first according to the following conditions Is each point a one-dimensional Morse critical point?

[0053] ;

[0054] Set a fixed-length statistical window along the path direction (window length is...). , (For the statistical window radius), calculate the local complexity density along the path. :

[0055] ;

[0056] in, This is an indicator function; it takes the value 1 when the condition within the parentheses is met, and 0 otherwise.

[0057] The local complexity density is assigned back to the corresponding coordinate point of the virtual conduction path in three-dimensional space to construct the distribution of topological disorder along the fiber conduction path. :

[0058] .

[0059] It should be noted that in an ideal cardiac physiological model, the propagation of the electrical excitation wave and the subsequent mechanical contraction wave along the myocardial fibers should be monotonically progressive, or contain only a very small number of modal changes (e.g., a single systolic-diastolic transition) throughout the entire systolic cycle. However, errors arising from multimodal fusion often manifest as high-frequency oscillations or local signal reversals along the fiber direction. This non-physical knotting phenomenon corresponds topologically to an abnormal increase in the number of scalar field critical points (i.e., maxima or minima). Therefore, by statistically analyzing the one-dimensional Morse critical point density along the virtual conduction path, a quantitative index can be constructed. This index does not depend on the absolute amplitude of the signal but purely reflects the structural complexity of the signal, thereby locating pathological or error regions that violate the laws of myocardial mechanical conduction.

[0060] Preferably, the calculation process first targets each virtual transmission path. One-dimensional activation function obtained by upsampling Differential analysis is performed. To obtain a stable derivative estimate on discrete data, the first derivative is calculated using the central difference method with second-order precision. The calculation formula is as follows: In this formula, Indicates the first The first path The rate of change of activation time along the fiber direction at each sampling point (i.e., the derivative). and These are the activation time values ​​for adjacent sampling points in the subsequent and preceding sequences, respectively. For a fixed step size for sampling along the path (e.g.) (millimeters). Next, Morse critical point detection is performed. This is done by employing the sign-flipping criterion, i.e., for any internal node... If the conditions are met Then determine the point. A one-dimensional Morse critical point exists nearby. Specifically, the sign of the derivative changes locally (from positive to negative or vice versa), corresponding to a peak or trough in the function curve. Subsequently, the local complexity density is calculated. A sliding statistical window is defined along the fiber direction, with a window radius of... (For example, take) Millimeters, a typical scale corresponding to the thickness of the myocardial wall. For each point on the path. Its local complexity density According to the formula Calculation. Among them, Indicates the arc length distance. The set of all neighboring points; Let be the characteristic function, when point The value is 1 when it is determined to be a critical point, and 0 otherwise; the denominator This represents the total length of the statistical window. Finally, the calculated density values... The value is then assigned back to the corresponding coordinates of the sampling point in three-dimensional space. This allows for the construction of a three-dimensional distribution of topological disorder along the fiber conduction pathway. If multiple paths pass through the same spatial voxel, then the arithmetic mean of the results calculated for these paths at that voxel is taken.

[0061] In some possible embodiments, a more smooth variant algorithm can be used when calculating the one-dimensional derivative and statistical density. For example, in cases with low signal-to-noise ratios, to avoid spurious critical points caused by digitized discrete noise, the one-dimensional sequence can be pre-processed before calculating the derivative. A five-point Savitzky-Golay filter is applied for preprocessing, or the smoothed derivative can be calculated directly using the Savitzky-Golay differential operator. Furthermore, when calculating the local complexity density, a Gaussian weighted window can be used instead of the rectangular window (BoxcarWindow) in the preferred embodiment described above. In this case, the density calculation formula is adjusted to... ,in The standard deviation of the Gaussian kernel (e.g., set to ) millimeters). The normalization coefficients are represented by this Gaussian weighting. This Gaussian weighting makes the calculated turbulence field more continuous and smooth in space, avoiding data jumps caused by window edge truncation effects, and is more conducive to subsequent use as continuous control parameters to guide adaptive filtering.

[0062] In a preferred embodiment, adjusting the filter parameters according to the distribution of topological disorder along the fiber conduction includes:

[0063] The distribution of topological disorder along the fiber conduction path is calculated according to the following formula. Linear normalization is performed to obtain the normalized disorder value. (Values ​​range from 0 to 1):

[0064] ;

[0065] in, Distribution of topological disorder along fiber conduction The global minimum value, Distribution of topological disorder along fiber conduction The global maximum value, To prevent small positive numbers from being divided by zero;

[0066] Construct a mapping function from the normalized disorder value to the local smoothness scale according to the following formula, and calculate the local smoothness scale at each point along the fiber direction. :

[0067] ;

[0068] in, This is the preset minimum smoothing scale along the fiber direction. The preset maximum smoothness scale along the fiber direction (satisfying) ), p is a preset nonlinear enhancement index (satisfying) This results in regions with higher topological disorder having larger local smoothness scales.

[0069] It is important to note that in cardiac biomechanical activation imaging, data quality varies significantly across different regions: in some regions, signals propagate smoothly along the fiber direction, and applying substantial smoothing filters can mask clinically significant local heterogeneity (such as conduction delays at scar margins); while in other regions, fusion errors or artifacts cause severe signal tangling along the fiber direction, requiring high-intensity smoothing to restore reasonable signal integrity. Therefore, by introducing a nonlinear mapping strategy, using topological disorder degree (FPIMC) as an inverse confidence indicator, the size of the filter kernel is dynamically modulated. The nonlinear enhancement term is used to construct a soft threshold gating system: for low-complexity regions (FPIMC close to 0), the smoothing scale is suppressed to near its minimum, achieving full-pass or weak noise reduction; as complexity increases, the smoothing scale rapidly and nonlinearly grows, strongly intervening in highly disordered regions. This ensures that the filtering operation only corrects non-physical structures, preserving the true physiological information to the greatest extent possible.

[0070] Preferably, the process of adjusting the filter parameters first involves performing global data normalization. This traverses the distribution of fiber conduction topological disorder across the entire left ventricular myocardium. Calculate its global minimum value and global maximum value Using the linear normalization formula Calculate the dimensionless normalized disorder value In this formula, It is a very small positive number (e.g.) ), only for use in equal To prevent numerical calculation errors such as zero denominators in extremely flat conditions and ensure program stability, the hyperparameters of the mapping function are determined. A minimum smoothing scale along the fiber direction is preset. The value is 1.0 mm, which typically corresponds to one voxel size or the smallest resolvable unit in an image, meaning that only discrete noise at the sub-voxel level is eliminated in the most reliable region; the maximum smoothing scale is preset. The value is 5.0 mm, set based on the average thickness of the left ventricular wall or the characteristic scale of the myocardial segment, meaning that the smoothing range can cover the entire local segment in the least reliable region; a pre-set nonlinear enhancement index is also included. for Finally, the local smoothing scale is calculated point by point. The calculation formula is: In this formula, The term utilizes the power function property to suppress the weight of the low-value range (e.g.) This makes the smooth scale of most low-disorder regions... Extremely close This protects the activation details of normal myocardium; and when Approaching At that time, this item rapidly increases, driving... Approaching This allows for powerful reshaping of highly disordered areas.

[0071] In some possible embodiments, other forms of monotonically increasing nonlinear functions can also be used when constructing the mapping function from disorder to smoothness. For example, a sigmoid function or a hyperbolic tangent function can be used instead of a power function. In a variant based on sigmoid, the mapping formula is designed as follows: ,in This is the steepness coefficient. This S-shaped mapping is... The most drastic change occurs at the intermediate value, while saturation occurs at both ends, providing a smooth control characteristic similar to a binary switch but maintaining continuous differentiability. Furthermore, regarding the range of normalization, in addition to global normalization, local adaptive normalization can also be used, i.e. and Take from the current point Statistical values ​​within a certain neighborhood range are used to accommodate the differences in basic complexity among different myocardial segments.

[0072] In a preferred embodiment, adaptive smoothing along only the virtual conduction path is performed on the whole left ventricular myocardial biomechanical activation distribution data to generate myocardial activation distribution data corrected for fiber structure constraints, including:

[0073] For each of the aforementioned virtual transmission paths Each sampling point on Based on the local smoothing scale at that point Construct a one-dimensional Gaussian weighted kernel that extends only along this path direction;

[0074] The one-dimensional Gaussian weighted kernel is used to calculate a weighted average of the full left ventricular myocardial biomechanical activation distribution data of the sampling point within the path neighborhood, resulting in a smoothed sampling point value. :

[0075] ;

[0076] in, Virtual transmission path Upper The numerical values ​​of the whole left ventricular myocardial biomechanical activation distribution data corresponding to each sampling point. The normalization coefficient is calculated using the following formula:

[0077] ;

[0078] Map the smoothed sampling point values ​​on all virtual transmission paths back to 3D spatial coordinates. If the same 3D spatial coordinate is covered by multiple paths, then take all the values ​​covering that coordinate. The average value is used to generate the myocardial activation distribution data corrected for fiber structure constraints. :

[0079] ;

[0080] in, To cover three-dimensional spatial coordinates The number of virtual transmission paths.

[0081] It should be noted that existing three-dimensional Gaussian filtering is typically performed in Cartesian coordinates. Applying smoothing equally along all three axes can lead to erroneous mixing across myocardial fiber layers, such as incorrectly contaminating late-activation regions of the epicardium with early activation signals from the endocardium. By locking the one-dimensional convolution onto the virtual conduction path, the filtering operation is effectively restricted to the propagation trajectory of the myocardial electromechanical waves. This approach precisely repairs wavefront breaks or non-physical oscillations propagating along the fiber direction while maintaining absolute sharpness in the direction perpendicular to the fiber, thus fully preserving the transmural gradients and interlaminar heterogeneity of the myocardium.

[0082] Preferably, this adaptive smoothing process is performed path-by-path based on the Lagrangian viewpoint. First, each virtual transmission path is traversed. (in And locate each discrete sampling point on the path. Read the local smoothing scale corresponding to this spatial point from the smoothing scale field. Next, using the current point... A one-dimensional discrete Gaussian weighted kernel is constructed centered on the path arc length and extending only along that direction. The corrected path numerical values... The result is obtained through the discrete convolution formula: In this calculation formula, Indicates the first The path is in the location The activation time value after smoothing correction; This indicates the path's location in the neighborhood. The original activation time observation at the location; and Both represent the arc length coordinates along the path; It represents the square of the geodesic distance between two points along the fiber trajectory; This is the local adaptive smoothing scale at that point. This determines the width of the Gaussian kernel smoothing range; This represents the step size for path sampling; For The cutoff window range centered on is usually taken as Coverage range. The denominator in the formula. It is the normalization coefficient that ensures energy conservation, and its calculation formula is: This coefficient ensures that the sum of the weights is 1, preventing the overall value from drifting due to smoothing operations. After completing the convolution calculations for all one-dimensional paths, it is necessary to process the data distributed along these paths. Write back to three-dimensional voxel space to generate the final field data. Since paths may intersect or overlap in space, for any voxel in a 3D mesh... If multiple paths pass through the voxel (or the voxel falls within the interpolation neighborhood of multiple path sampling points), then collect the correction values ​​on all relevant paths. The arithmetic mean of the values ​​is calculated as the final myocardial activation distribution data for that voxel, corrected for fiber structure constraints. .

[0083] In some possible implementations, other forms of kernel functions or write-back strategies can be used when performing path smoothing. For example, in addition to Gaussian kernels, one-dimensional bilateral filter kernels can be used. In this variation, the weights depend not only on the spatial distance. It also depends on the intensity difference of the activation time value. The formula is expanded to This bilateral filtering strategy smooths out clutter while better protecting genuine wavefront transitions (such as conduction blockage at the edge of scar tissue), preventing excessive blurring. Furthermore, during the data writing-back (Splatting) stage, if multiple paths overlap, in addition to simple arithmetic averaging, a confidence-weighted average can be used. That is, for each point on a path, a weight is assigned based on the signal-to-noise ratio of its original data or the degree of matching between the path's tangent vector at that point and the local fiber direction. Paths with higher matching degrees or higher signal-to-noise ratios contribute more weight. Value in the final The larger the proportion, the better.

[0084] In a preferred embodiment, calculating the corrected residual conduction disorder distribution includes:

[0085] Based on the myocardial activation distribution data corrected by fiber structure constraints As the input object, it is again constructed along the virtual conduction path built by the three-dimensional myocardial fiber architecture model. Sampling yields the corrected one-dimensional activation function. ;

[0086] Calculation using the finite difference method The derivative with respect to arc length The critical point is determined by the derivative sign detection, and the corrected local complexity density is obtained statistically. :

[0087] ;

[0088] in, Let be the radius of the complexity statistical window along the fiber direction. This is an indicator function; it takes the value 1 when the condition within the parentheses is met, and 0 otherwise.

[0089] Will Mapping back to the corresponding coordinates in three-dimensional space generates the residual conduction disorder distribution. This distribution is used to quantify the degree of non-physical topological back-reversal along the fiber direction that still exists after adaptive smoothing:

[0090] .

[0091] It should be noted that although the full left ventricular myocardial biomechanical activation distribution data has been structurally smoothed based on the initial topological disorder, in complex clinical pathology (such as severe post-infarction scar junctions) or under extreme imaging artifact conditions, simple smoothing may not completely eliminate all non-physical reentry, or the smoothed data may still retain topological uncertainties in some high-gradient regions. Therefore, by modifying the corrected data... A Morse analysis was performed again along the fiber path to quantify the degree of remaining persistent topological knots. The resulting distribution of residual conduction disorder was then analyzed. It is no longer just a filter control parameter, but has been transformed into a terminal diagnostic indicator that reflects the integrity of the data structure or clinical reliability.

[0092] Preferably, the process of calculating the residual conduction disorder distribution directly uses the geometric information of the virtual conduction path to ensure the consistency of the analytical benchmark and reduce computational overhead. First, myocardial activation distribution data corrected for fiber structure constraints is used. For the input object, along the first Virtual transmission path Resampling was performed. Because... The function has already been updated in three-dimensional space; the sampled result is the corrected one-dimensional activation function. Subsequently, the derivative calculation and critical point detection process is performed on the new function. Specifically, the corrected path derivative is calculated. And detect sign reversal. For any point on the path. If satisfied If the condition is met, then the point is determined to be a corrected one-dimensional Morse critical point. Next, the corrected local complexity density is calculated using a sliding window. Its calculation formula is In this formula, The statistical window radius is typically kept consistent with the value used when calculating the initial disorder (e.g., (millimeters) to facilitate comparison before and after; Represents the set of neighborhood indexes covered by the window; This is a critical point indicator function for the corrected data. Finally, the calculated density values... Mapping back to 3D space coordinates Generate a three-dimensional residual conduction disorder distribution. In this process, if a certain region, after undergoing strong smoothing, has a function curve along the fiber direction... If it becomes monotonous or contains only a single extreme value that conforms to physiological laws, then that place If the value approaches zero, it indicates that the data structure has returned to normal; conversely, if... If the value remains high, it indicates that the data in that area contains deep structural contradictions that cannot be repaired through smoothing, and thus belongs to a low-confidence area.

[0093] In some possible embodiments, the relative improvement rate can be introduced as an auxiliary indicator when calculating the residual disorder, or a higher-order analysis based on curvature can be used. In a variant embodiment based on the relative improvement rate, the calculation of the residual distribution considers not only the corrected absolute complexity but also the initial complexity. Apply weighting. Define the relative residual index. This indicator can highlight persistent noise regions that should have been smoothed out but still exist. In another variation, considering that the smoothed curve may already be very smooth and no longer exhibit obvious sign reversal (i.e., the first derivative disappears at zero-crossing points), but there may still be non-physiological sharp curvature (i.e., abnormal second derivative), the critical point detection can be replaced with or expanded to inflection point detection, that is, statistical second derivative. The zero-crossing density serves as a more sensitive indicator of residual structure anomalies. Furthermore, to improve computational efficiency, in scenarios where extremely high spatial resolution is not required, resampling can be performed only on the initial disorder level, without resampling all paths. High-risk paths exceeding a certain threshold are locally recalculated, while the residual disorder in other areas is set to zero by default.

[0094] In a preferred embodiment, the residual conduction disorder distribution is mapped to the myocardial activation distribution data corrected by fiber structure constraints in a multidimensional visual manner to generate a heart failure follow-up visualization image that simultaneously presents the activation sequence and its structural reliability on different visual channels, including:

[0095] Establish mapping function This enables the calculation of three-dimensional points from the left ventricular geometric model. To the polar coordinate system of the two-dimensional bullseye diagram The mapping relationship projects data in three-dimensional space onto a two-dimensional pixel plane;

[0096] Calculate each pixel Corresponding three-dimensional region The average value within, where This represents the projected mean of myocardial activation distribution data after correction for fiber structure constraints. The projected mean of the residual conduction disorder distribution:

[0097] ;

[0098] in, Represents a three-dimensional region The number of three-dimensional points contained within;

[0099] To each and Performing linear normalization yields values ​​ranging from 0 to 1. (Normalized average of myocardial activation distribution data corrected for fiber structure constraints) and (Normalized average distribution of residual conduction disorder);

[0100] Each pixel is encoded using the Hue-Saturation-Luminosity (HSV) color model, with the following formula:

[0101] ;

[0102] ;

[0103] ;

[0104] in, For the hue components of a pixel, by Linearity determines the timing of activation. and This is a preset hue range; For the saturation component of a pixel, by The negative exponential function determines that regions with higher residual disorder have lower saturation, thus representing a decrease in structural reliability. This is the preset saturation decay coefficient; For the luminance component of a pixel, This is a preset fixed brightness constant.

[0105] It's important to note that in existing Bull's-eye charts or 3D renderings, color typically encodes only a single parameter (such as activation time). Doctors cannot determine whether the displayed values ​​are based on high-quality signals or are in a state of high noise or topological disorder. By introducing the HSV (Hue-Saturation-Luminosity) color model, single parameters are mapped to hue, and structural reliability (the inverse of residual disorder) is mapped to saturation. Thus, in areas with smooth data structures that conform to the physiological conduction patterns of the myocardium, the image displays vibrant colors, facilitating doctor readings; while in areas with chaotic data structures and severe non-physical backlash, the image becomes grayish-white, visually reducing the salience of unreliable information and providing a direct warning and preventing misleading interpretations.

[0106] Preferably, the multidimensional visual mapping process first involves a geometric projection from a three-dimensional anatomical space to a two-dimensional parametric space. This is achieved by establishing a geometric model of the left ventricle. To the polar coordinate system of the two-dimensional bullseye diagram One-to-one mapping function Typically, this mapping places the apex of the left ventricle at the origin of polar coordinates, the base at the circumference of the largest radius circle, and the anterior, septal, inferior, and lateral walls evenly distributed at angles. To generate high-resolution pixel-level images, a reverse ray projection method is employed: for each discrete pixel on the two-dimensional bullseye plane... Determine its corresponding region in three-dimensional space under the inverse transformation of the mapping function. (Typically, it's a ray passing through the myocardial wall or a thin cone). Next, the raw channel values ​​for that pixel are calculated. The corresponding region is then statistically analyzed. Myocardial activation distribution data with fiber structure constraint correction for all three-dimensional voxels contained herein. and residual conduction disorder distribution The arithmetic mean of the values ​​is calculated as the activation time value of that pixel. and complexity value Then, numerical normalization is performed. For activation time, the formula is used. Map it to The interval, in which and The statistical extreme values ​​of the entire graph; for complexity, use the formula Map it to The range is then determined. Finally, HSV encoding calculation is performed. Hue components are then considered. The activation time is linearly determined, and the calculation formula is as follows: Specifically, the settings (Blue indicates the earliest activation) and (Red, representing the latest activation), making the color change conform to the intuitive understanding of warm and cool tones. Saturation component. The residual disorder is determined by a negative exponential function, and the calculation formula is as follows: In this formula, The preset saturation attenuation coefficient is preferably set to [value]. . Used to control the sensitivity of visual cues: when When the value is large, even a small residual disturbance can cause a sharp drop in saturation, and the image will quickly turn gray. Luminance component Set as a fixed constant (For example or This ensures that the overall image has sufficient brightness. Through this encoding, the color of each pixel in the final generated image is determined. Both incorporate information about time and quality.

[0107] In some possible implementations, this can be applied to discretized visualization scenarios based on the AHA17 segment model, or other color mixing spaces can be used. In a variant of the AHA17 segment model, instead of calculating for each pixel, the 3D data is first aggregated into the average of 17 standard segments, and then these 17 discrete blocks are encoded using the HSV encoding described above. The image generated in this way consists of 17 sectors of different colors. Although the spatial resolution is reduced, it is more in line with the format of standard clinical reports. In the color space variant, transparency (AlphaChannel) can be used instead of saturation as a carrier of credibility. In this case, the background layer is preset to a gray or black grid, and the transparency of the foreground layer is... From the formula Calculation. When the disorder level is high, the foreground layer becomes transparent, revealing the gray background underneath, achieving the same visual effect of desaturating or blurring unreliable areas. Furthermore, the above HSV encoding logic also applies to 3D surface rendering; simply change the computational object from two-dimensional pixels. Replace with 3D mesh vertices This allows the generation of a 3D heart model with a textured surface that conveys structural credibility.

[0108] like Figure 2 As shown, Figure 2 It demonstrates how to construct an analytical foundation that conforms to the anatomical characteristics of the heart and interprets the geometric relationship between the three-dimensional myocardial fiber architecture model and the virtual conduction path. Figure 2The array of blue arrows distributed within the left ventricular wall represents a three-dimensional model of the myocardial fiber architecture, following the continuous spiral change of myocardial fibers from the endocardium (positive helix angle) to the adventitia (negative helix angle). The red, yellow, green, and other multi-colored curves interspersed within represent virtual conduction paths. These paths are not arbitrary lines, but rather streamlines generated by numerical integration strictly along the blue fiber vector field, starting from seed points at specific anatomical locations. Figure 2 The complex and anisotropic three-dimensional cardiac body data is transformed into a set of one-dimensional signals along fiber tracks that conform to the physiological electromechanical wave conduction laws.

[0109] like Figure 3 As shown, Figure 3 The complete algorithmic flow of topological disorder along the fiber conduction path, from the original signal to the quantified index, is explained in detail. The top subplot shows a one-dimensional activation function (i.e., the curve of activation time versus path arc length) sampled on a specific virtual conduction path, which exhibits non-monotonic fluctuations in certain regions. The middle subplot calculates the spatial derivative of this function and marks the one-dimensional Morse critical points (i.e., local peaks or troughs) where the derivative sign flips. These points physically correspond to anomalous signal foldbacks or non-physical oscillations. The bottom subplot shows the local complexity density, a stepped distribution generated by statistically analyzing the density of the aforementioned critical points within a sliding window along the path. Figure 3 By capturing extreme points that should not exist along the fiber direction, signal noise that is difficult to visually measure is transformed into a quantitative score of topological disorder.

[0110] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A method for visualizing and processing heart failure follow-up data based on multimodal image fusion, characterized in that: include: A geometric model of the left ventricle and a three-dimensional myocardial fiber architecture model were constructed based on multimodal follow-up images, and initial full left ventricular myocardial biomechanical activation distribution data were generated by fusing them. A virtual conduction path is constructed along the three-dimensional myocardial fiber architecture model, and the distribution of fiber conduction topological disorder, characterizing the degree of non-physical re-entry of the activation wave along the fiber, is calculated, including: For each one-dimensional activation function on the virtual conduction path, the derivative of its arc length is calculated using the finite difference method. The sign change of the derivative is detected. When the product of the derivatives of the preceding and following adjacent points is negative, the point is determined to be a one-dimensional Morse critical point. A fixed-length statistical window is set along the path direction, and the total number of one-dimensional Morse critical points contained in the statistical window is counted and divided by the window length to obtain the local complexity density. The local complexity density is used as the value of the corresponding coordinate point of the virtual conduction path in three-dimensional space to construct the distribution of topological disorder along the fiber conduction. Based on the distribution of fiber conduction topology disorder, the filtering parameters are adjusted, and adaptive smoothing is performed only along the virtual conduction path on the whole left ventricular myocardial biomechanical activation distribution data to generate myocardial activation distribution data corrected by fiber structure constraints. The corrected residual conduction disorder distribution is calculated and mapped in a multidimensional visual manner with the myocardial activation distribution data corrected by fiber structure constraints, generating a heart failure follow-up visualization image that simultaneously presents the activation sequence and structural reliability on different visual channels.

2. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 1, characterized in that, A left ventricular geometric model and a three-dimensional myocardial fiber architecture model were constructed based on multimodal follow-up images, including: A reference geometric time point is selected in the multimodal follow-up images, and a left ventricular reference geometric grid is established through image segmentation. The normalized transmural depth at each point in the left ventricular reference geometric grid is calculated, where the normalized transmural depth is the linearly normalized distance from the endocardium to the epicardium. Based on a preset helical fiber model, the endocardial helical angle and the epicardial helical angle are set, and the helical angle is linearly interpolated according to the normalized transmural depth to obtain the local helical angle at each point. Based on the local cylindrical coordinate system, the circumferential direction vector is weighted by the cosine value of the local helical angle, and the vertical direction vector is weighted by the sine value of the local helical angle, and the three-dimensional myocardial fiber architecture model is synthesized and normalized.

3. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 2, characterized in that, Initial full-ventricular myocardial biomechanical activation distribution data were generated by fusing the left ventricular geometric model and the three-dimensional myocardial fiber architecture model, including: Mechanical activation time data from multiple different imaging modalities are registered onto the left ventricular geometric model to form a single-modal activation distribution field corresponding to each modality. A fixed fusion weight is pre-set for each imaging modality, wherein the fixed fusion weight is a non-negative number and the sum is one. Weighted linear summation is performed on the single-modal activation distribution fields corresponding to each modality to calculate the whole left ventricular myocardial mechanical activation distribution data.

4. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 3, characterized in that, Constructing a virtual conduction path along the three-dimensional myocardial fiber architecture model includes: Several representative seed points are selected within the left ventricular geometric model, and these seed points are evenly distributed within the left ventricular myocardial wall. Starting from each seed point, the model extends bidirectionally along the direction indicated by the three-dimensional myocardial fiber architecture model using a numerical integration method with a fixed step size until it exceeds the myocardial boundary, generating a virtual conduction path in three-dimensional space. Along each virtual conduction path, the whole left ventricular myocardial biomechanical activation distribution data is sampled at fixed intervals, and the three-dimensional data is transformed into a one-dimensional activation function distributed along the arc length of the path.

5. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 4, characterized in that, Adjusting the filter parameters according to the aforementioned distribution of topological disorder along the fiber conduction path includes: The minimum and maximum values ​​of the fiber conduction topological disorder distribution within the entire left ventricle are obtained. The fiber conduction topological disorder distribution is then linearly normalized to obtain the normalized disorder. A minimum smoothing scale and a maximum smoothing scale along the fiber direction, as well as a nonlinear enhancement exponent greater than one, are pre-defined. A mapping function from the normalized disorder to the local smoothing scale is constructed. This mapping function linearly maps the normalized disorder value to the minimum and maximum smoothing scales after exponentially operating the nonlinear enhancement exponent, thereby making the local smoothing scale larger for regions with higher topological disorder.

6. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 5, characterized in that, The whole left ventricular myocardial biomechanical activation distribution data is subjected to adaptive smoothing along the virtual conduction path only to generate myocardial activation distribution data corrected by fiber structure constraints, including: For each sampling point on each virtual conduction path, a one-dimensional Gaussian weighted kernel extending along the direction of the virtual conduction path is constructed based on the local smoothing scale at that sampling point. Using the one-dimensional Gaussian weighted kernel, a weighted average is calculated on the whole left ventricular myocardial biomechanical activation distribution data in the neighborhood of that sampling point on the virtual conduction path to obtain the smoothed sampling point value. The smoothed sampling point values ​​on all virtual conduction paths are mapped back to three-dimensional spatial coordinates to generate the myocardial activation distribution data corrected by fiber structure constraints.

7. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 6, characterized in that, The calculation of the corrected residual conduction disorder distribution includes: Using the myocardial activation distribution data corrected by fiber structure constraints as input, a virtual conduction path is constructed again along the three-dimensional myocardial fiber architecture model, and the corrected one-dimensional activation function is extracted. The derivative sign detection and critical point statistics are performed again on the corrected one-dimensional activation function to calculate the corrected local complexity density. The corrected local complexity density is mapped back to three-dimensional space to generate the residual conduction disorder distribution. The residual conduction disorder distribution is used to quantitatively characterize the degree of non-physical backtracking of the topological structure along the fiber direction that still exists after adaptive smoothing.

8. The method for visualizing heart failure follow-up data based on multimodal image fusion according to claim 7, characterized in that, The residual conduction disorder distribution is mapped in a multidimensional visual manner to the myocardial activation distribution data corrected by fiber structure constraints, generating a heart failure follow-up visualization image that simultaneously presents the activation sequence and structural reliability on different visual channels, including: A mapping relationship is established from the left ventricular geometric model to the two-dimensional bullseye polar coordinate system, projecting the data in three-dimensional space onto the two-dimensional pixel plane. For each pixel on the two-dimensional pixel plane, the average value of the myocardial activation distribution data corrected by fiber structure constraints and the average value of the residual conduction disorder distribution in the corresponding three-dimensional region are calculated. The average values ​​are linearly normalized respectively. Each pixel is encoded using a hue-saturation-luminance color model: the hue component of the pixel is linearly determined by the normalized average value of the activation distribution data to express the timing of activation; the saturation component of the pixel is determined by the normalized average value of the residual disorder distribution through a negative exponential function, so that the region with higher residual disorder has lower saturation to express the decrease in structural reliability; the luminance component of the pixel is set to a fixed constant.

Images