Fetal echocardiography system for identifying and assessing valvular motion abnormalities

Abstract

This invention relates to a system for identifying and assessing valvular motion abnormalities in fetal echocardiography, specifically in the field of echocardiography imaging. The system first robustly extracts the valvular central skeleton from a single frame of ultrasound image, overcoming image noise interference. Then, it constructs a physically constrained curve model to simulate the valvular tissue's mechanical properties and precisely tracks its motion trajectory in the image sequence through spatiotemporal optimization. Finally, the system isolates the overall cardiac motion, quantifies the valve's own opening and closing deformation, interleaflet motion synchronicity, and local tissue strain rate, thereby achieving early, objective, and quantitative accurate assessment of fetal valvular motion dysfunction and assisting in clinical diagnosis.

Description

Technical Field

[0001] This invention relates to the field of cardiac ultrasound imaging, and more specifically, to a system for identifying and assessing valvular motion abnormalities in fetal cardiac ultrasound. Background Technology

[0002] In prenatal screening and diagnosis of congenital heart disease in fetuses, the structural and functional assessment of heart valves is crucial. Valve morphological abnormalities or functional insufficiency, such as stenosis or regurgitation, are key indicators or accompanying lesions of many complex cardiac malformations. In clinical practice, two-dimensional grayscale echocardiography is the primary imaging method for assessing the fetal heart. Doctors observe the opening and closing of valves through dynamic images. However, current assessments largely rely on the physician's subjective visual judgment and experience, lacking objective and precise quantitative analysis of valve motion patterns. Valvees, especially the leaflets of the aortic and mitral valves, are not... Rigid bodies undergo rapid and complex non-rigid motions during each cardiac cycle, including bending, stretching, and twisting. Early or mild valvular dysfunction often manifests first as a reduction in local motion amplitude, subtle asynchrony between different leaflets, or abnormal motion at specific locations (such as the valve tip compared to the valve root). However, in ultrasound sequences with significant image noise, limited frame rate, and interference from the overall fetal motion, it is difficult for the human eye to stably and accurately capture and measure these subtle spatiotemporal changes, leading to insufficient sensitivity in identifying early lesions and potentially delaying diagnosis.

[0003] Currently, although some studies have attempted to quantify valve motion using image processing techniques, existing methods have significant limitations. Traditional optical flow methods rely on the assumption of constant image brightness, which often fails in ultrasound images due to speckle noise and variations in the acoustic reflection characteristics of tissue interfaces. This leads to computational failures in valve regions with indistinct texture features, and motion estimation results are easily contaminated by noise. Block matching methods are ill-suited for non-rigid deformations and have limited accuracy. Both optical flow and block matching output motion fields that are superficial and hybrid motions, including the relative motion of the valve itself with respect to the heart, as well as the translation of the heart as a whole within the thoracic cavity. Regarding the relationship between valve opening and closing and rotation, existing technologies lack effective decoupling mechanisms and cannot isolate the "effective motion" component that purely reflects the valve's opening and closing function. More importantly, these methods typically treat the valve as a whole or track a few points, failing to characterize its continuous deformation field as a flexible membrane. Consequently, they cannot quantify key mechanical parameters such as local strain, curvature changes, and time delays in motion between different spatial points. How to robustly estimate the high-precision motion trajectory of each point on the valve surface from low-quality, two-dimensional ultrasound image sequences, and separate the quantitative indicators of local deformation and motion desynchronization from them, is a core problem that existing technologies have not yet solved. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a system for identifying and assessing valvular motion abnormalities in fetal cardiac ultrasound. The system utilizes an initialization module, a modeling module, an optimization module, and an analysis module to solve the problems mentioned in the background section.

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: it includes an initialization module, a modeling module, an optimization module, and an analysis module connected in sequence, wherein;

[0006] Initialization module: Used to receive ultrasound image sequence of fetal heart, select a frame image from the ultrasound image sequence, perform morphological opening operation on a set of linear structuring elements with increasing scale for the single frame image, calculate the morphological gradient of the image after opening operation at each scale, preserve the ridge structure that is stable at multiple consecutive scales through scale space analysis, and generate and output a skeleton with a single pixel width by performing a final erosion process that preserves the connection of the union of the preserved stable ridge structures.

[0007] Modeling module: It receives a skeleton with a width of one pixel, detects the branch points and endpoints of the skeleton to use them as fixed key control points, and adaptively inserts additional control points on the skeleton segments between the key points according to the local curvature changes, thereby generating a set of control points. Based on this set of control points, it constructs a physical elastic curve model containing spring damper connections and discrete bending energy terms to complete the initialization of the physical elastic curve model.

[0008] The optimization module receives the ultrasound image sequence and the initialized physical elastic curve model. It constructs a global time-domain energy function and minimizes the function using a multi-resolution optimization strategy to obtain the control point motion trajectory of the physical elastic curve model over the entire image sequence time range. The global time-domain energy function includes data terms, internal constraint terms, and time-domain smoothing terms.

[0009] Analysis module: It is used to receive the motion trajectory of the control point, decompose the effective motion displacement relative to the cardiac stability reference point from the motion trajectory of the control point, and calculate the curvature weighted time delay between the leaflets based on the time when the effective motion of each control point reaches the peak and the local curvature of its motion trajectory to quantify motion desynchronization. At the same time, it obtains the motion vector field of the valve surface by interpolation based on the effective motion displacement of all control points and calculates the local strain rate.

[0010] In a preferred embodiment, the initialization module performs morphological opening operations sequentially using a set of scale-increasing linear structuring elements, and calculates the morphological gradient of the image after the opening operations at each scale. Specifically, this includes:

[0011] Select a frame from the received ultrasound image sequence and define a set of scale-increasing linear structuring elements;

[0012] For each structural element at each scale, a morphological opening operation is first performed on the selected image to obtain the opening operation image at that scale.

[0013] Subsequently, for each opening operation image, a first intermediate image is obtained by performing a dilation operation using a preset compact structuring element with a size smaller than all linear structuring elements, and a second intermediate image is obtained by performing an erosion operation. The values ​​of corresponding pixels in the first intermediate image and the second intermediate image are subtracted to obtain the morphological gradient image corresponding to that scale, thereby generating a set of morphological gradient images that correspond one-to-one with multiple scales.

[0014] The ridge structure that is stable across multiple consecutive scales is preserved through scale-space analysis. The specific operation is as follows:

[0015] For a set of morphological gradient images, for each pixel position in a frame of an image, the following judgment and calculation are performed: determine the scale at which the gradient value of the pixel position first exceeds a preset gradient threshold in a set of morphological gradient images, and take it as the starting scale; starting from the starting scale, along the scale increment direction, count the total number of adjacent scales in which the gradient value of the pixel position continuously exceeds the gradient threshold, and define the total number as the scale persistence measure of the pixel position.

[0016] Simultaneously, a set of gradient values ​​is obtained at these continuous scales where the gradient value at the pixel location continuously exceeds the gradient threshold. The maximum and minimum values ​​in this set of gradient values ​​are found, and the difference between the maximum and minimum values ​​is calculated to obtain the range of the set of gradient values. Then, the ratio of the range to a corrected denominator is calculated. The corrected denominator is the sum of the maximum value in the set of gradient values ​​and a preset minimum positive constant. Finally, the ratio is subtracted from the number 1, and the result is defined as the gradient consistency measure at the pixel location.

[0017] A minimum continuous scale number threshold and a gradient consistency threshold are preset. For each pixel position in a frame of image, the pixel is determined to belong to a ridge structure that is stable at multiple consecutive scales only if the scale persistence metric of the pixel position is not less than the minimum continuous scale number threshold and the gradient consistency metric of the pixel position is higher than the gradient consistency threshold.

[0018] All pixel locations identified as belonging to stable ridge structures are marked to generate a binary stable ridge marker image.

[0019] In a preferred embodiment, after generating the stable ridge marker image, the union of the preserved stable ridge structures is subjected to a final erosion process that preserves connectivity. The specific process is as follows:

[0020] First, morphological closing operations are performed on the stable ridge marker image to connect the small gaps in the image caused by noise or insufficient contrast, resulting in a aggregated ridge image.

[0021] Calculate the morphological skeleton of the aggregated ridge image to obtain a preliminary skeleton image composed of single-pixel wide foreground pixels;

[0022] Calculate the distance transformation map of the aggregated ridge image. The value of each foreground pixel in the distance transformation map represents the distance from that point to the nearest background pixel in the image.

[0023] Detect skeleton branch point pixels in the preliminary skeleton image. Skeleton branch point pixels refer to pixels in the preliminary skeleton image that are adjacent to three or more other skeleton pixels.

[0024] Based on the distance transformation map and the positions of the detected skeleton branch point pixels, a connectivity-preserving iterative morphological erosion operation is performed on the aggregated ridge image. The rules of the iterative morphological erosion operation are as follows: In each iteration, only foreground pixels in the aggregated ridge image that simultaneously meet the following two conditions are allowed to be removed: First, the pixel is not located within a protected neighborhood centered on any skeleton branch point pixel with a preset pixel distance as its radius; Second, the distance value corresponding to the pixel in the distance transformation map is greater than the number of erosion iterations already performed.

[0025] Repeatedly perform the iterative morphological erosion operation until the width of all remaining foreground connected regions in the aggregated ridge image becomes one pixel, or until a preset maximum number of iterations is reached. At this point, stop the iteration and use the resulting image as a single-pixel-wide skeleton.

[0026] In a preferred embodiment, after obtaining the single-pixel-wide skeleton, a skeleton post-processing operation is further included before output, which specifically includes:

[0027] Branch analysis and cleanup are performed on the single-pixel wide skeleton image obtained after the final erosion process that preserves connectivity.

[0028] The branch analysis and cleanup are achieved through the following steps: identifying all endpoint pixels in the single-pixel wide skeleton image, where an endpoint pixel is a pixel in the single-pixel wide skeleton image that is only adjacent to one other skeleton pixel.

[0029] For each identified endpoint pixel, starting from that endpoint pixel, traverse along the path of the skeleton pixels connected to it, and record all the pixels traversed until a skeleton branch point pixel or another endpoint pixel is encountered. The skeleton branch point pixel refers to a pixel in the single-pixel wide skeleton image that is connected to three or more other skeleton pixels.

[0030] The total number of pixels traversed from the starting endpoint pixel to the ending point is defined as the length of the skeleton branch originating from that starting endpoint pixel.

[0031] A preset length threshold is set; if the calculated length of a certain skeleton branch is less than the length threshold, the skeleton branch is determined to be a short branch caused by noise, and all pixels constituting the skeleton branch are deleted from the single-pixel wide skeleton image.

[0032] After performing the above traversal, length calculation, judgment and pixel deletion operations on all identified endpoint pixels, the cleaned skeleton containing only the main trunk and significant branches with a width of one pixel is finally obtained and output.

[0033] In a preferred embodiment, the modeling module detects the branch points and endpoints of the skeleton as fixed key control points, and adaptively inserts additional control points on the skeleton segments between the key points based on local curvature changes. The specific process is as follows:

[0034] First, traverse the received single-pixel-wide skeleton, identify all skeleton branch point pixels and endpoint pixels, and define all identified skeleton branch point pixels and endpoint pixels together as a set of fixed key control points.

[0035] Subsequently, on the skeleton path segment determined by any two adjacent fixed key control points, additional control points are adaptively inserted based on the local curvature values ​​of each pixel on the path segment. The calculation process of the local curvature value is as follows: For a pixel to be calculated on the path segment, obtain the position coordinates of its direct predecessor pixel and direct successor pixel on the path segment, calculate the first vector from the predecessor pixel to the pixel to be calculated, and the second vector from the pixel to be calculated to the successor pixel, and calculate the angle between the first vector and the second vector; then calculate the sum of the magnitudes of the first vector and the second vector; finally, multiply the value by the absolute value of the sine of the angle, and divide by the sum to obtain the local curvature value of the pixel to be calculated.

[0036] After calculating the local curvature values ​​of each pixel on the skeleton path segment, the local curvature values ​​are compared with preset high curvature thresholds and medium curvature thresholds. The comparison process is as follows: in path segments where the local curvature value is greater than the high curvature threshold, additional control points are inserted at a first preset interval; in path segments where the local curvature value is between the medium curvature threshold and the high curvature threshold, additional control points are inserted at a second preset interval greater than the first preset interval; in path segments where the local curvature value is less than the medium curvature threshold, additional control points are inserted at a third preset interval greater than the second preset interval, or no additional control points are inserted.

[0037] The additional control points generated through the above comparison process constitute an additional control point set. The union of the fixed key control point set and the additional control point set is the generated set of control points. The initial position of each control point is set to the coordinates of the corresponding pixel in the single-pixel-wide skeleton.

[0038] In a preferred embodiment, the specific process of constructing a physical elastic curve model based on the set of control points, including spring-damper connections and discrete bending energy terms, is as follows:

[0039] Each point in a set of control points is defined as a mass point, and adjacent points on a single-pixel-wide skeleton in this set of control points are connected using spring-damper units to form a chain structure. After the connection is completed, based on the geometric information of the single-pixel-wide skeleton, the parameters in the physical elastic curve model are initialized. The initialization includes:

[0040] Q1. For each pair of adjacent control points connected by spring damper units, define a natural length parameter, the value of which is set to the shortest path pixel distance between the two control points in a single-pixel-wide skeleton multiplied by a preset pixel physical size.

[0041] Q2. For each pair of adjacent control points connected by spring damper units, define a spring stiffness coefficient, which is set to be proportional to the reciprocal of the local or average width of the skeleton in the corresponding region.

[0042] Q3. Define a bending stiffness coefficient for any three consecutive control points in a chain structure. For control points located near the skeleton branch point pixels in the fixed key control point set, the associated bending stiffness coefficient is set to a higher value than that of the non-branch point region.

[0043] Q4. For any three control points arranged consecutively in a chain structure, define a balance angle parameter. Its value is directly set to the exterior angle value calculated by the three control points based on their initial position coordinates in a single-pixel-wide skeleton.

[0044] Q5. For each control point or each spring damper connection, define a damping coefficient that is proportional to the product of the square root of the corresponding spring stiffness coefficient and the mass parameter assigned to the control point.

[0045] After completing the above initialization assignment, define the internal elastic potential energy of the physical elastic curve model. The internal elastic potential energy consists of two parts: spring tension energy and discrete bending energy.

[0046] By defining the internal elastic potential energy and initializing the parameters as described above, the physical elastic curve model is initialized.

[0047] In a preferred embodiment, the global time-domain energy function constructed in the optimization module is specifically constructed as follows:

[0048] First, the optimization module receives the ultrasound image sequence and the initialized physical elastic curve model. The state of the physical elastic curve model at any time is defined by the position of all its control points at that time. The solution objective is to obtain the motion trajectory of all control points of the physical elastic curve model over the entire time range of the ultrasound image sequence.

[0049] To solve for the motion trajectory, a global time-domain energy function is constructed. This function is a weighted sum of the data terms, internal constraint terms, and time-domain smoothing terms. Specifically, the global time-domain energy function equals the sum of the following three terms:

[0050] The first step is to calculate the data items for each moment in the ultrasound image sequence from the first moment to the last moment, and then sum all the data items.

[0051] The second term involves multiplying an internal constraint term weight coefficient by the internal constraint term for each moment in the ultrasound image sequence from the first moment to the last moment, calculating the internal constraint term for that moment, and summing all these internal constraint terms.

[0052] The third term is a time-domain smoothing term weighting coefficient multiplied by each time point from the second time point to the last time point in the ultrasound image sequence, calculating the time-domain smoothing term between that time point and the previous time point, and summing all these time-domain smoothing terms.

[0053] The calculation process for the data items is as follows:

[0054] First, for each frame of the ultrasound image sequence, calculate its phase consistency feature map and local principal direction field;

[0055] Then, for the continuous curve defined by all control point positions at a certain moment in the physical elasticity curve model, a series of points are sampled on the curve; for each sampled point, the following operations are performed:

[0056] The feature intensity value at the sampling point is extracted from the phase consistency feature map. Based on the local shape of the continuous curve near the sampling point, the tangent direction of the curve at the point is calculated and normalized into a vector of unit length, called the unit tangent vector of the sampling point. The local principal direction angle value of the image structure at the sampling point is extracted from the local principal direction field and converted into a unit direction vector. The absolute value of the angle between the unit tangent vector and the unit direction vector of the sampling point is calculated. A Gaussian kernel function with a preset standard deviation parameter is applied to the absolute value to map it, and a direction consistency weight value between zero and one is obtained.

[0057] Finally, the feature intensity value is multiplied by the orientation consistency weight value to obtain the orientation matching contribution value of the sampling point; the orientation matching contribution values ​​of all sampling points are summed and negatively taken as the main part of the data item at that moment. This data item also includes a distance penalty sub-item. The final value of the data item is obtained by adding the sum of the negative orientation matching contribution values ​​to the value of the distance penalty sub-item.

[0058] The calculation process of the internal constraint term is as follows: at each moment, the internal elastic potential energy of the physical elastic curve model at that moment is directly calculated. The internal elastic potential energy is composed of the spring tension energy and the discrete bending energy. The value of the internal constraint term is equal to the sum of the spring tension energy and the discrete bending energy.

[0059] The calculation process of the time-domain smoothing term is as follows: For any two adjacent time points, denoted as the previous time point and the current time point, the time-domain smoothing term is composed of a weighted sum of two parts: a smoothness penalty term based on motion acceleration and a penalty term based on motion pattern consistency.

[0060] The calculation process for the smoothness penalty term based on motion acceleration is as follows: For each control point in the physical elastic curve model, firstly, calculate its displacement vector from the previous moment to the current moment. This displacement vector is equal to the control point's position coordinate at the current moment minus its position coordinate at the previous moment. Simultaneously, calculate the control point's displacement vector relative to the even earlier moment from the previous moment. Then, calculate the difference between the two displacement vectors to obtain an approximate acceleration vector of the control point at the current moment. Next, calculate the square of the magnitude of this approximate acceleration vector, multiply it by a first acceleration weighting coefficient, and obtain the first acceleration weighting coefficient of the control point. First, an acceleration penalty value is calculated. Simultaneously, the ratio of the approximate magnitude of the acceleration vector to a correction denominator is calculated, and this ratio is squared and multiplied by a second acceleration weighting coefficient to obtain the second acceleration penalty value for that control point. The correction denominator is the sum of the magnitude of the control point's displacement vector at the previous moment and 0.001. Finally, a preset control point weighting coefficient, the first acceleration penalty value, and the second acceleration penalty value are multiplied together to obtain the acceleration smoothing penalty contribution value for that control point. The acceleration smoothing penalty contribution values ​​for all control points are summed to obtain the smoothness penalty term based on motion acceleration.

[0061] The calculation process for the motion pattern consistency-based penalty term is as follows: First, the sets of displacement vectors of all control points at the previous time step and the sets of displacement vectors of all control points at the current time step are statistically analyzed. Each set is considered as an empirical probability distribution of a displacement vector. Then, the Kullback-Leibler divergence of the empirical distribution of the displacement vector at the current time step relative to the empirical distribution of the displacement vector at the previous time step is calculated. Finally, the Kullback-Leibler divergence is multiplied by a pattern consistency weight coefficient to obtain the motion pattern consistency-based penalty term.

[0062] The value of the time-domain smoothing term is equal to the sum of the smoothness penalty term based on motion acceleration and the penalty term based on motion pattern consistency.

[0063] In a preferred embodiment, the multi-resolution optimization strategy specifically includes:

[0064] First, based on the initial static positions of each control point in the initialized physical elastic curve model received from the modeling module, a preliminary control point position sequence that is uniformly distributed in time and runs through the entire ultrasound image sequence is generated using a linear interpolation method, and this preliminary sequence is defined as the initial motion trajectory guess.

[0065] Next, Gaussian pyramid downsampling is performed on the ultrasound image sequence and the initial motion trajectory guess in the time dimension. Specifically, low-pass filtering and interval sampling are performed on both along the time axis to construct a multi-layer time resolution data from coarse to fine, including the downsampled ultrasound image sequence and the downsampled control point state sequence.

[0066] Then, an iterative optimization process is executed, starting from the coarsest time resolution layer and progressing layer by layer to the finest layer, specifically:

[0067] R1. At the current time resolution layer, based on the downsampled ultrasound image sequence corresponding to this layer and the control point state sequence of the current layer, construct and minimize the constructed global time domain energy function, and solve to obtain the motion trajectory of the control point of the current layer.

[0068] R2. If the current layer is not the finest original temporal resolution layer, the control point motion trajectory obtained in step R1 is upsampled to the next layer with a finer temporal resolution through interpolation on the time axis, and this upsampled trajectory is used as the initial control point state sequence of the next layer optimization process.

[0069] R3. Set the next layer as the new current layer and repeat R1 and R2.

[0070] The optimization process starts from the coarsest temporal resolution layer and iterates through steps R1 to R3 until the global temporal energy function is minimized at the finest original temporal resolution layer. Finally, the physical elastic curve model is output as the control point motion trajectory over the entire time range of the original ultrasound image sequence.

[0071] In a preferred embodiment, the analysis module decomposes the effective motion displacement from the control point trajectory and calculates the curvature-weighted time delay between the leaflets. The specific process is as follows:

[0072] First, based on the received motion trajectory of the control point, the motion trajectory of a cardiac stability reference point is determined; the cardiac stability reference point is obtained by automatically tracking a preset, relatively fixed anatomical landmark in the ultrasound image sequence; then, the motion trajectory of the cardiac stability reference point is subjected to time low-pass filtering to extract its overall motion component that changes slowly over time.

[0073] Based on this overall motion component, for each moment and each control point in the physical elastic curve model, a motion decomposition operation is performed:

[0074] S1. Calculate the absolute position change of the control point from the initial time to the current time. The change is the position coordinate of the control point at the current time minus its position coordinate at the initial time.

[0075] S2. Calculate a rigid motion displacement that the control point should undergo, induced by the overall motion of the cardiac stability reference point. This rigid motion displacement is obtained through the following steps:

[0076] S2.1 Obtain the relative position vector of the control point relative to the initial overall motion component of the cardiac stability reference point at the initial moment. The vector is the position coordinate of the control point at the initial moment minus the coordinate of the overall motion component of the reference point at the initial moment.

[0077] S2.2 Obtain an optimal rotation matrix, which is obtained by solving a minimization problem. The goal of the minimization problem is to minimize the sum of squares of the total difference between the position of all control points relative to the current overall motion component of the reference point at the current time and the result obtained by applying the optimal rotation matrix to the relative position vector of each control point relative to the initial overall motion component of the reference point at the initial time.

[0078] S2.3 Apply the optimal rotation matrix to the relative position vector obtained in step S2.1;

[0079] S2.4 Add the result of step S2.3 to the displacement vector of the overall motion component of the reference point from the initial time to the current time. The displacement vector is the overall motion component coordinate of the reference point at the current time minus its overall motion component coordinate at the initial time.

[0080] S2.5. Subtract the position coordinates of the control point at the initial moment from the result obtained in step S2.4 to finally obtain the rigid motion displacement;

[0081] S3. Subtract the rigid motion displacement calculated in step S2 from the absolute position change calculated in step S1. The result is the effective motion displacement of the control point at that moment.

[0082] Through the above processes S1 to S3, the effective motion displacement of all control points at all times is obtained, forming an effective displacement vector field.

[0083] After obtaining the effective displacement vector field, the curvature-weighted time delay is calculated.

[0084] In a preferred embodiment, the specific process of obtaining the motion vector field of the valve surface and calculating the local strain rate based on the effective motion displacement interpolation of all control points is as follows:

[0085] First, a continuous motion vector field is constructed based on the effective displacement vector field. The process is as follows: For each moment in the ultrasound image sequence, the static position coordinates of all control points of the physical elastic curve model at the initial moment are used as the source point for interpolation, and the effective motion displacement vectors corresponding to each control point at the same moment are used as the target value of the source point. A two-dimensional vector field is fitted on all source points using the thin-plate spline interpolation method. Specifically, this is achieved by minimizing an objective function consisting of a data fitting term and a regularization term. The data fitting term is calculated as the sum of the squares of the Euclidean norms of the differences between the two-dimensional vector field function values ​​to be fitted at all source points and the target effective displacement vectors corresponding to those points. The regularization term is used to penalize excessive curvature of the vector field function. Its value is a preset regularization parameter multiplied by the result of an integral expression calculated over the entire two-dimensional image plane. By solving the problem of minimizing the objective function, a continuous motion vector function defined on the two-dimensional image plane is obtained. The value of this function at any point represents the effective motion velocity vector of that point at the current moment, thus obtaining the continuous motion vector field of the valve surface.

[0086] After obtaining the continuous motion vector field, the local strain rate is calculated as follows: For each location point within the continuous motion vector field, the following steps are performed sequentially:

[0087] T1. Calculate the spatial velocity gradient tensor at this location: Decompose the velocity of the continuous motion vector field at the current point into horizontal and vertical components; calculate the partial derivatives of the horizontal velocity component with respect to the horizontal coordinate and the vertical coordinate, respectively; these four partial derivative values ​​together form a 2x2 matrix, called the spatial velocity gradient tensor at this point; the partial derivatives are approximated on a regular grid using the central difference method;

[0088] T2. Calculate the strain rate tensor at this location: Add the spatial velocity gradient tensor obtained in step T1 to its own transpose matrix, and then multiply the result by one-half. The result is called the strain rate tensor at this point. This tensor is a symmetric second-order tensor, and its components represent the local tensile or compressive rates of the medium at this point in different directions.

[0089] T3. Extract one or more scalarized local strain rate indices from the strain rate tensor, specifically including:

[0090] T3.1 Maximum Principal Strain Rate: The eigenvalues ​​of the strain rate tensor obtained in step T2 are used as the maximum principal strain rate at that point. This value characterizes the rate of maximum elongation or shortening in the local deformation at that point.

[0091] T3.2 Equivalent tensile and compressive strain rate: The sum of the squares of the components of the strain rate tensor obtained in step T2 is multiplied by two-thirds, and then the square root of the product is taken as the equivalent tensile and compressive strain rate at that point. This value comprehensively reflects the overall deformation rate intensity caused by pure tension or compression at that point.

[0092] T3.3 Shear strain rate: Calculate the determinant value of the strain rate tensor obtained in step T2, and calculate the square of its trace; when the trace of the strain rate tensor is not zero, divide four times the determinant value by the square of the trace, and then take the square root of the quotient value. The result is taken as the shear strain rate at that point.

[0093] By sequentially performing steps T1 to T3 on all locations on the valve surface, the maximum principal strain rate distribution map, the equivalent tensile-compressive strain rate distribution map, and the shear strain rate distribution map covering the entire valve region are obtained. These distribution maps together constitute the local strain rate.

[0094] The beneficial effects of this invention are as follows: First, the scheme robustly extracts the valve centerline skeleton from a single frame of ultrasound image, overcoming image noise interference; then, it constructs a physically constrained curve model to simulate the mechanical properties of valve tissue, and accurately tracks its motion trajectory in the image sequence through spatiotemporal optimization; finally, the scheme isolates the overall motion of the heart, quantifies the valve's own opening and closing deformation, the synchronicity of movement between valve leaflets, and the strain rate of local tissues, thereby achieving an early, objective, and quantitative accurate assessment of fetal heart valve motion dysfunction, and assisting in clinical diagnosis. Attached Figure Description

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

[0096] Figure 2 This is a block diagram of the system structure of the present invention. Detailed Implementation

[0097] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0098] In the description of this application, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the stated features. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0099] In the description of this application, the term "for example" is used to mean "used as an example, illustration, or description." Any embodiment described as "for example" in this application is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is provided to enable any person skilled in the art to make and use the invention. Details are set forth in the following description for purposes of explanation. It should be understood that those skilled in the art will recognize that the invention can be made without using these specific details. In other instances, well-known structures and processes will not be described in detail to avoid obscuring the description of the invention with unnecessary detail. Therefore, the invention is not intended to be limited to the embodiments shown, but is consistent with the broadest scope of the principles and features disclosed in this application.

[0100] Example 1

[0101] This embodiment provides, for example Figure 1-2The system shown is a fetal cardiac ultrasound system for identifying and assessing valvular motion abnormalities, specifically comprising: an initialization module, a modeling module, an optimization module, and an analysis module connected in sequence;

[0102] Initialization module: Used to receive ultrasound image sequence of fetal heart, select a frame image from the ultrasound image sequence, perform morphological opening operation on a set of linear structuring elements with increasing scale for the single frame image, calculate the morphological gradient of the image after opening operation at each scale, preserve the ridge structure that is stable at multiple consecutive scales through scale space analysis, and generate and output a skeleton with a single pixel width by performing a final erosion process that preserves the connection of the union of the preserved stable ridge structures.

[0103] Modeling module: It receives a skeleton with a width of one pixel, detects the branch points and endpoints of the skeleton to use them as fixed key control points, and adaptively inserts additional control points on the skeleton segments between the key points according to the local curvature changes, thereby generating a set of control points. Based on this set of control points, it constructs a physical elastic curve model containing spring damper connections and discrete bending energy terms to complete the initialization of the physical elastic curve model.

[0104] The optimization module receives the ultrasound image sequence and the initialized physical elastic curve model. It constructs a global time-domain energy function and minimizes the function using a multi-resolution optimization strategy to obtain the control point motion trajectory of the physical elastic curve model over the entire image sequence time range. The global time-domain energy function includes data terms, internal constraint terms, and time-domain smoothing terms.

[0105] Analysis module: Receives the motion trajectories of control points, decomposes the effective motion displacement relative to the cardiac stability reference point from these trajectories, and calculates the curvature-weighted time delay between valve leaflets based on the peak time of the effective motion at each control point and the local curvature of its trajectory to quantify motion desynchronization. Simultaneously, it interpolates the motion vector field on the valve surface based on the effective motion displacements of all control points and calculates the local strain rate.

[0106] In this embodiment, it is specifically noted that in the initialization module, a set of linear structuring elements with increasing scale are used to perform morphological opening operations sequentially, and the morphological gradient of the image after the opening operation at each scale is calculated. This specifically includes:

[0107] Select a frame from the received ultrasound image sequence and define a set of linear structuring elements with increasing scale, where the length of each structuring element is set according to the typical width range of the fetal valve (for example, the length can be set to 3, 5, 7, 9 pixels in sequence to cover the thinnest to the thickest possible valve structure).

[0108] For each structural element at each scale, a morphological opening operation is first performed on the selected image to obtain the opening operation image at that scale (this operation can effectively filter out noise points and small artifacts that are narrower than the current structural element, while preserving linear structures such as valves).

[0109] Subsequently, for each opening operation image, a first intermediate image is obtained by dilation operation using a pre-defined compact structuring element (such as a 3×3 pixel square) smaller than all linear structuring elements, and a second intermediate image is obtained by erosion operation. The values ​​of corresponding pixels in the first intermediate image and the second intermediate image are subtracted to obtain the corresponding morphological gradient image at that scale, thereby generating a set of morphological gradient images that correspond one-to-one with multiple scales (the morphological gradient images highlight significant region boundaries at that observation scale).

[0110] The ridge structure that is stable across multiple consecutive scales is preserved through scale-space analysis. The specific operation is as follows:

[0111] For a set of morphological gradient images, for each pixel position in a frame of an image, perform the following judgment and calculation:

[0112] The starting scale is determined when the gradient value of a pixel in a set of morphological gradient images first exceeds a preset gradient threshold (the gradient threshold can be set according to the overall contrast and noise level of the image, for example, 10%-20% of the maximum gray value of the image). Starting from the starting scale, along the scale increment direction, the total number of adjacent scales in which the gradient value of the pixel continuously exceeds the gradient threshold is counted, and this total number is defined as the scale persistence measure of the pixel.

[0113] Simultaneously, a set of gradient values ​​is obtained at these continuous scales where the gradient value at the pixel position continuously exceeds the gradient threshold. The maximum and minimum values ​​in this set of gradient values ​​are found, and the difference between the maximum and minimum values ​​is calculated to obtain the range of this set of gradient values. Then, the ratio of the range to a correction denominator is calculated. The correction denominator is the sum of the maximum value in this set of gradient values ​​and a preset minimum positive constant (the minimum positive constant is used to prevent the denominator from being zero, and can usually be taken as 1e-6).

[0114] Finally, subtract the ratio from the number 1 and define the result as the gradient consistency measure of the pixel location (the closer the measure is to 1, the smoother the gradient response of the point at each scale, and the more likely it is to belong to the real structure rather than noise fluctuations).

[0115] A minimum continuous scale number threshold and a gradient consistency threshold are preset (for example, the minimum continuous scale number threshold can be set to 3 and the gradient consistency threshold can be set to 0.7). For each pixel location in a frame of image, the pixel is determined to belong to a ridge structure that is stable across multiple continuous scales only if the scale persistence metric of the pixel location is not less than the minimum continuous scale number threshold and the gradient consistency metric of the pixel location is higher than the gradient consistency threshold (this dual criterion can effectively distinguish between real, stable anatomical edges and random noise or transient artifacts).

[0116] All pixel locations identified as belonging to stable ridge structures are labeled to generate a binary stable ridge label image, wherein the labeled pixel locations have a first pixel value (e.g., 1), and the unlabeled pixel locations have a second pixel value (e.g., 0) that is different from the first pixel value.

[0117] After generating a stable ridge marker image, the union of the preserved stable ridge structures is subjected to a final erosion process that preserves connectivity. The specific process is as follows:

[0118] First, a morphological closing operation is performed on the stable ridge marker image (using 3x3 or 5x5 square or circular structuring elements) to connect the small gaps in the image caused by noise or insufficient contrast, resulting in a aggregated ridge image (this design ensures that the valve contours, which should be continuous, are topologically connected).

[0119] Calculate the morphological skeleton of the aggregated ridge image to obtain a preliminary skeleton image composed of single-pixel wide foreground pixels;

[0120] Calculate the distance transformation map of the aggregated ridge image. The value of each foreground pixel in the distance transformation map represents the distance from that point to the nearest background pixel in the image (this map provides spatial information for subsequent conditional erosion).

[0121] The process involves detecting skeleton branch point pixels in the preliminary skeleton image. Skeleton branch point pixels are pixels that are adjacent to three or more other skeleton pixels in the preliminary skeleton image (these points typically correspond to critical topological structures such as valve junctions and require special protection). Based on the distance transform map and the positions of the detected skeleton branch point pixels, a connectivity-preserving iterative morphological erosion operation is performed on the aggregated ridge image. The rules for the iterative morphological erosion operation are as follows:

[0122] In each iteration, only foreground pixels in the aggregated ridge image that meet both of the following conditions are allowed to be removed: First, the pixel is not located within the protected neighborhood centered on any skeleton branch point pixel and with a preset pixel distance as the radius (the preset pixel distance, i.e. the protected radius, can usually be set to 2 pixels).

[0123] Second, the distance value of the pixel in the distance transformation map is greater than the number of erosion iterations that have been performed so far (this condition ensures that erosion starts from the outermost layer and proceeds layer by layer inward, and pixels with a distance value greater than the number of iterations are located in the current "layer to be eroded").

[0124] Repeatedly perform iterative morphological erosion operations until the width of all remaining foreground connected regions in the aggregated ridge image becomes one pixel, or a preset maximum number of iterations is reached (e.g., 10 times to prevent infinite looping in abnormal situations). At this point, stop the iteration and use the resulting image as a single-pixel-wide skeleton (this method effectively avoids the topological connection breakage that may be caused by traditional ultimate erosion while achieving skeleton refinement, especially protecting the structural integrity at branch points).

[0125] After obtaining the skeleton with a width of one pixel, post-processing of the skeleton is performed before output, which specifically includes:

[0126] Branch analysis and cleanup are performed on the single-pixel wide skeleton image obtained after the final erosion process that preserves connectivity.

[0127] The specific process of branch analysis and cleanup is as follows: Identify all endpoint pixels in the single-pixel wide skeleton image. An endpoint pixel is a pixel in the single-pixel wide skeleton image that is only adjacent to one other skeleton pixel.

[0128] For each identified endpoint pixel, starting from that endpoint pixel, traverse along the path of the skeleton pixels connected to it, and record all the pixels traversed until a skeleton branch point pixel or another endpoint pixel is encountered. A skeleton branch point pixel refers to a detected pixel that connects to three or more other skeleton pixels in a single-pixel wide skeleton image.

[0129] The total number of pixels traversed from the starting endpoint pixel to the ending point (skeleton branch point pixel or another endpoint pixel) is defined as the length of the skeleton branch originating from that starting endpoint pixel.

[0130] A preset length threshold is set (this length threshold can be estimated and set according to the image resolution and the physical size of the valve, for example, 5 pixels in length; branches shorter than this length are considered meaningless noise). If the length of the calculated skeleton branch is less than the length threshold, the skeleton branch is determined to be a short branch caused by noise, and all pixels constituting the skeleton branch are deleted from the single-pixel wide skeleton image.

[0131] After performing the above traversal, length calculation, judgment and pixel deletion operations on all identified endpoint pixels, the cleaned skeleton containing only the main trunk and significant branches is finally obtained and output (this step removes burrs and isolated noise lines that may remain from the previous processing, ensuring the purity of the output skeleton and providing high-quality topology input for the subsequent modeling module).

[0132] In this embodiment, it is specifically noted that in the modeling module, branch points and endpoints of the skeleton are detected and used as fixed key control points. Additional control points are then adaptively inserted into the skeleton segments between these key points based on local curvature changes. The specific process is as follows:

[0133] First, traverse the received single-pixel-wide skeleton, identify all skeleton branch point pixels and endpoint pixels, and define all identified skeleton branch point pixels and endpoint pixels together as a set of fixed key control points.

[0134] Subsequently, on the skeleton path segment determined by any two adjacent fixed key control points, additional control points are adaptively inserted based on the local curvature values ​​at each pixel point on the path segment; the calculation process for the local curvature values ​​is as follows:

[0135] For a pixel to be calculated on a path segment, obtain the position coordinates of its immediate predecessor and immediate successor pixels on that path segment. Calculate the first vector from the predecessor pixel to the pixel to be calculated, and the second vector from the pixel to the successor pixel, and calculate the angle between the first and second vectors. Then calculate the sum of the magnitudes of the first and second vectors. Finally, multiply the value by the absolute value of the sine of the angle, and divide by the sum. The result is the local curvature value of the pixel to be calculated. This calculation method, by simultaneously considering the turning angle of the path point and the adjacent step size, can more accurately characterize the minute bending changes of the skeleton at the pixel level. Compared with the method of simply calculating the angle between adjacent line segments, it is less sensitive to noise and pixel discretization effects.

[0136] After calculating the local curvature values ​​of each pixel on the skeleton path segment, the local curvature values ​​are compared with preset high curvature thresholds and medium curvature thresholds. The comparison process is as follows:

[0137] In path segments where the local curvature value is greater than the high curvature threshold, additional control points are inserted at a first preset interval.

[0138] In path segments where the local curvature value is between the medium curvature threshold and the high curvature threshold, additional control points are inserted at a second preset interval distance greater than the first preset interval distance. By segmenting according to curvature value and inserting additional control points at different densities, higher geometric description accuracy can be provided in key areas with high valve curvature (such as the valve tip) without significantly increasing the total number of control points in the model (i.e., computational complexity). This allows for better capture of subtle deformations in these areas during subsequent motion tracking. In this application, the high curvature threshold can be set to 0.8, the medium curvature threshold can be set to 0.5, and the specific values ​​of the first preset interval distance, the second preset interval distance, and the third preset interval distance can be adjusted according to the image resolution and the desired model accuracy.

[0139] In path segments where the local curvature value is less than the medium curvature threshold, insert additional control points at a third preset interval distance that is greater than the second preset interval distance, or do not insert additional control points.

[0140] The additional control points generated through the above comparison process constitute an additional control point set. The union of the fixed key control point set and the additional control point set is the generated set of control points. The initial position of each control point is set to the coordinates of the corresponding pixel in a single-pixel-wide skeleton.

[0141] Based on this set of control points, a physical elastic curve model including spring-damper connections and discrete bending energy terms is constructed. The specific process is as follows:

[0142] Each point in a set of control points is defined as a mass point, and adjacent points on a single-pixel-wide skeleton in this set of control points are connected using spring-damper units to form a chain structure. After the connection is completed, based on the geometric information of the single-pixel-wide skeleton, the parameters in the physical elastic curve model are initialized. The initialization includes:

[0143] Q1. For each pair of adjacent control points connected by spring-damper units, define a natural length parameter. Its value is set to the shortest path pixel distance between the two control points in a single-pixel-wide skeleton multiplied by a preset pixel physical size. This definition enables the initial equilibrium state of the model to accurately fit the actual meandering path of the single-pixel-wide skeleton, rather than a simple straight-line distance, thereby injecting the correct geometric information into the model at the start of the motion.

[0144] Q2. For each pair of adjacent control points connected by spring damper units, define a spring stiffness coefficient. This coefficient is set to be proportional to the reciprocal of the local or average width of the skeleton in the corresponding region. This setting is based on a simplified biomechanical prior: regions with thicker tissue generally have higher stiffness. By relating the spring stiffness coefficient to the local width, the model can roughly simulate the stiffness differences that may exist in different parts of the valve (such as the valve body and valve tip).

[0145] Q3. Define a bending stiffness coefficient for any three consecutive control points in the chain structure. For control points located near the skeleton branch point pixels in the fixed set of key control points, the associated bending stiffness coefficient is set to a higher value than that of non-branch point regions. At the valve junction (corresponding to the skeleton branch point), the tissue structure is usually more stable and has relatively low mobility. By appropriately increasing the bending stiffness coefficient at this point, the relative rigidity of this junction can be simulated in the physical model, preventing unreasonable excessive bending in subsequent optimization.

[0146] Q4. For any three control points arranged consecutively in the chain structure, define a balance angle parameter. Its value is directly set to the exterior angle value calculated based on the initial position coordinates of the three control points in the single-pixel-wide skeleton. This balance angle parameter defines the model's "memory shape", that is, the angle to which the model tends to recover when no external force is applied. This parameter is calculated directly using the geometry of the initial skeleton, ensuring that the model's internal balance state is consistent with the observed initial anatomical morphology of the valve.

[0147] Q5. For each control point or each spring-damper connection, define a damping coefficient that is proportional to the square root of the corresponding spring stiffness coefficient and the mass parameter assigned to the control point. This proportionality aims to provide the system with a coefficient close to the critical damping, which helps to suppress unreasonable oscillations in subsequent dynamic simulations or optimizations, making the model's motion response smoother, more stable, and more consistent with the viscoelastic characteristics of biological tissues.

[0148] After completing the above initialization assignment, define the internal elastic potential energy of the physical elastic curve model. The internal elastic potential energy consists of two parts: spring tension energy and discrete bending energy.

[0149] The calculation process for spring tension energy is as follows: For each pair of adjacent control points directly connected by spring damper units in the chain structure, perform the following operations:

[0150] A1. Calculate the current Euclidean distance between the control points;

[0151] A2. Subtract the natural length parameter defined in Q1 from the current Euclidean distance to obtain the elongation;

[0152] A3. Calculate the square of the elongation;

[0153] A4. Multiply the squared value by half of the spring stiffness coefficient defined by Q2 to obtain the basic tensile energy of the spring unit;

[0154] A5. Calculate a correction factor, the value of which is: one plus an anharmonicity coefficient multiplied by the square of the elongation, and then divided by the sum of the square of the natural length parameter defined in Q1 and a preset minimum positive constant; this correction factor is introduced to simulate the "strain hardening" effect exhibited by biological soft tissue under large tension, that is, the stiffness increases with the increase of tension, which makes the mechanical behavior of the model closer to the real valve tissue characteristics than a simple linear spring (Hooke's Law);

[0155] A6. Multiply the basic tensile energy by the correction factor to obtain the final spring tensile energy contributed by the spring damper unit.

[0156] Repeat A1 to A6 for all spring damper units, and add up all the final spring tension energies obtained. The sum is the spring tension energy.

[0157] The calculation process for discrete bending energy is as follows: For any three control points arranged consecutively in a chain structure, perform the following operations:

[0158] B1. Obtain the position coordinates of these three control points;

[0159] B2. Based on the position coordinates of these three points, calculate the vector from the middle point to the point before it, and the vector from the middle point to the point after it.

[0160] B3. Calculate the dot product of the two vectors above, and calculate the magnitude of each vector respectively;

[0161] B4. Divide the dot product by the product of the two moduli and calculate the arccosine of the result to obtain an exterior angle value;

[0162] B5. Calculate the absolute value of the difference between the exterior angle value and the balance angle parameter defined in Q4 to obtain the deviation.

[0163] B6. Calculate the square of the deviation and multiply it by half of the bending stiffness coefficient defined by Q3 to obtain the bending energy contributed by these three points.

[0164] Repeat B1 to B6 for all possible consecutive three control points, and add up all the bending energies obtained. The sum is the discrete bending energy.

[0165] The sum of the spring's tensile energy and discrete bending energy is the internal elastic potential energy of the physical elastic curve model. By defining the internal elastic potential energy and initializing the parameters as described above, the physical elastic curve model is initialized. Thus, a valve dynamics model with complete physical properties (mass, elasticity, bending resistance, damping) and an initial equilibrium state is established. This model not only encodes the valve's geometry but, more importantly, imposes mechanical behavior constraints on it through physical parameters. This ensures that when solving its motion trajectory in subsequent image sequences, the results not only satisfy image matching but also conform to biomechanical rationality. This is the key to achieving high-precision and robust motion estimation.

[0166] In this embodiment, it is specifically necessary to explain the global time-domain energy function constructed in the optimization module, and the specific process is as follows:

[0167] First, the optimization module receives the ultrasound image sequence and the initialized physical elastic curve model. The state of the physical elastic curve model at any time is defined by the position of all its control points at that time. The solution objective is to obtain the motion trajectory of all control points of the physical elastic curve model over the entire time range of the ultrasound image sequence.

[0168] To solve for the motion trajectory, a global time-domain energy function is constructed. This function is a weighted sum of the data terms, internal constraint terms, and time-domain smoothing terms. Specifically, the global time-domain energy function equals the sum of the following three terms:

[0169] The first step is to calculate the data items for each moment in the ultrasound image sequence from the first moment to the last moment, and then sum all the data items.

[0170] The second term involves multiplying an internal constraint term weight coefficient by the internal constraint term for each moment in the ultrasound image sequence from the first moment to the last moment, calculating the internal constraint term for that moment, and summing all these internal constraint terms.

[0171] The third term, a weighted coefficient of a time-domain smoothing term, is multiplied by the time-domain smoothing term between each time point from the second time point to the last time point in the ultrasound image sequence, and the sum of all these time-domain smoothing terms is calculated. Among them, the data term is used to measure the degree of matching between the shape of the physical elastic curve model at any time point and the ultrasound image at the corresponding time point; the internal constraint term is used to constrain the deformation of the physical elastic curve model to conform to its preset physical elastic properties; the time-domain smoothing term is used to constrain the temporal continuity and smoothness of the motion of all control points in the physical elastic curve model.

[0172] The calculation process for the data items is as follows:

[0173] First, for each frame of the ultrasound image sequence, its phase consistency feature map and local principal orientation field are calculated. The phase consistency feature map reflects the stable position of structural features in the image and is insensitive to changes in illumination and contrast. The local principal orientation field describes the angle of the most significant structural orientation at each pixel in the image. In addition, the phase consistency feature map is calculated by processing the image at multiple scales and directions using a Log-Gabor filter bank, which can effectively suppress common speckle noise in ultrasound images and highlight stable structural features such as tissue edges and textures. The local principal orientation field can be obtained by calculating the principal eigenvector direction of the image gradient structure tensor.

[0174] Then, for the continuous curve defined by all control point positions at a certain moment in the physical elasticity curve model, a series of points are sampled on the curve; for each sampled point, the following operations are performed:

[0175] The feature intensity value at the sampling point is extracted from the phase consistency feature map. Based on the local shape of the continuous curve near the sampling point, the tangent direction of the curve at that point is calculated and normalized into a vector of unit length, called the unit tangent vector of the sampling point. The local principal direction angle value of the image structure at the sampling point is extracted from the local principal direction field and converted into a unit direction vector. The absolute value of the angle between the unit tangent vector and the unit direction vector of the sampling point is calculated. A Gaussian kernel function with a preset standard deviation parameter is applied to this absolute value to obtain a direction consistency weight value between zero and one. The Gaussian kernel function maps zero angle to a maximum value of one, and the weight value decreases as the angle increases. The preset standard deviation parameter controls the tolerance of direction consistency. The larger the value, the lighter the penalty for angle deviation. The preset standard deviation parameter can be taken in the range of 0.1 to 5.0. For example, it can be set to 1.0 to achieve a good balance between direction sensitivity and robustness.

[0176] Finally, the feature intensity value is multiplied by the orientation consistency weight value to obtain the orientation matching contribution value of the sampling point. The orientation matching contribution values ​​of all sampling points are summed and negative, which is used as the main part of the data item at that moment. The negative operation makes the process of minimizing the energy function equivalent to maximizing the consistency between the model curve and the image features in position and orientation. Optionally, the data item also includes a distance penalty sub-item. The calculation process of the sub-item is as follows: First, using the same or similar image processing methods as those used in the initialization module to generate stable ridge marker images and subsequent processing, an image representing the potential centerline position of the valve is extracted from the ultrasound image at the corresponding moment, which is defined as a multi-scale stable ridge map. Then, for each point sampled on the physical elastic curve model, the nearest distance from that point to the multi-scale stable ridge map is calculated. The Euclidean distance of the foreground pixels, or the distance value corresponding to the location of the sampling point in the distance transformation map of the multi-scale stable ridge map, is used as the spatial distance between the sampling point and the multi-scale stable ridge map. Finally, the spatial distances of all sampling points are summed and multiplied by a preset distance penalty weight coefficient to obtain the value of the distance penalty sub-item. The final value of the data item is obtained by adding the sum of the negative values ​​of the direction matching contribution to the value of the distance penalty sub-item. The distance penalty sub-item, as a weak global position prior, can prevent the model from drifting excessively due to lack of image force in frames with extremely weak or missing image texture information, and ensure the stability of the optimization process throughout the sequence. The distance penalty weight coefficient is usually set to a small value, such as 0.1, to avoid excessive interference with the main data item based on phase consistency and orientation field.

[0177] The calculation process for the internal constraint terms is as follows: At each moment, the internal elastic potential energy of the physical elastic curve model at that moment is directly calculated. This internal elastic potential energy consists of the sum of the spring tension energy and the discrete bending energy. The calculation method for the spring tension energy is exactly the same as that defined in the modeling module. That is, for each pair of adjacent control points connected by spring dampers in the model, the elongation of the current distance relative to the natural length is calculated, and the calculation is performed based on this elongation, the spring stiffness coefficient, and the anharmonicity correction factor. The calculation method for the discrete bending energy is also exactly the same as that defined in the modeling module. That is, for each pair of adjacent control points connected by spring dampers in the model, the elongation of the current distance relative to the natural length is calculated, and the calculation is performed based on this elongation, the spring stiffness coefficient, and the anharmonicity correction factor. For three consecutive control points, calculate the deviation of their current exterior angle from the equilibrium angle, and perform calculations based on this deviation and the bending stiffness coefficient; the value of the internal constraint term is equal to the sum of the spring tension energy and the discrete bending energy; the introduction of the internal constraint term ensures that the optimized motion trajectory not only fits the image, but also conforms to the physical deformation law of the valve tissue as an elastic membrane (such as resisting excessive stretching and bending), thus obtaining a more biomechanically reliable motion estimate. The weight coefficient of the internal constraint term is used to balance image matching and physical rationality, and its value needs to be adjusted experimentally, for example, it can be initialized to 1.0;

[0178] The calculation process of the time-domain smoothing term is as follows: For any two adjacent time points, denoted as the previous time point and the current time point, the time-domain smoothing term is composed of a weighted sum of two parts: a smoothness penalty term based on motion acceleration and a penalty term based on motion pattern consistency.

[0179] The smoothness penalty term based on motion acceleration is calculated as follows: For each control point in the physical elastic curve model, firstly, its displacement vector from the previous moment to the current moment is calculated. This displacement vector is equal to the control point's position coordinate at the current moment minus its position coordinate at the previous moment. Simultaneously, the displacement vector of the control point relative to the previous moment is calculated. Then, the difference between the two displacement vectors is calculated to obtain an approximate acceleration vector of the control point at the current moment. Next, the square of the magnitude of the approximate acceleration vector is calculated and multiplied by a first acceleration weighting coefficient to obtain the first acceleration penalty value of the control point. Simultaneously, the ratio of the magnitude of the approximate acceleration vector to a correction denominator is calculated, and the square of this ratio is multiplied by a second acceleration weighting coefficient to obtain the second acceleration penalty value of the control point. The correction denominator is the sum of the magnitude of the displacement vector of the control point at the previous moment and 0.001. This 0.001 is used to prevent division by zero or numerical instability in the relative acceleration calculation when the magnitude of the displacement vector at the previous moment is zero (or close to zero). Finally, a preset control point weight coefficient, the first acceleration penalty value, and the second acceleration penalty value are multiplied to obtain the acceleration smoothing penalty contribution value of the control point. The acceleration smoothing penalty contribution values ​​of all control points are summed to obtain the smoothness penalty term based on motion acceleration. By simultaneously constraining the absolute acceleration and the rate of change of relative acceleration, high-frequency jitter of the trajectory can be effectively suppressed, and the physiological characteristics of cardiac motion velocity changes (such as rapid systole and gradual diastole) can be adapted. The first acceleration weight coefficient can be set to 0.3, and the second acceleration weight coefficient can be set to 0.2.

[0180] The calculation process for the motion pattern consistency-based penalty term is as follows: First, the sets of displacement vectors of all control points at the previous time step and the sets of displacement vectors of all control points at the current time step are statistically analyzed. Each set is considered as an empirical probability distribution of displacement vectors. Then, the Kullback-Leibler divergence of the empirical distribution of displacement vectors at the current time step relative to the empirical distribution of displacement vectors at the previous time step is calculated. This Kullback-Leibler divergence measures the difference between the two probability distributions; the smaller the value, the more similar the two distributions are. Finally, the Kullback-Leibler divergence is multiplied by a pattern consistency weight coefficient to obtain the motion pattern consistency-based penalty term. This penalty term encourages the overall motion pattern of the valve (i.e., the statistical distribution of motion velocities of all parts) to remain stable between adjacent frames, avoiding unreasonable, global motion pattern abrupt changes, and enhancing the coherence and physiological rationality of long-term motion trajectories.

[0181] The value of the temporal smoothing term is equal to the sum of the smoothness penalty term based on motion acceleration and the penalty term based on motion pattern consistency. The weight coefficient of the temporal smoothing term controls the smoothness of the motion trajectory in the time dimension. If its value is too large, it will cause motion estimation to lag, and if it is too small, it will not be able to effectively suppress the jitter caused by noise. It needs to be adjusted according to the image frame rate.

[0182] The multi-resolution optimization strategy employed specifically includes:

[0183] First, based on the initial static positions of each control point in the initialized physical elastic curve model received from the modeling module, a preliminary control point position sequence that is uniformly distributed in time and runs through the entire ultrasound image sequence is generated using a linear interpolation method. This preliminary sequence is defined as the initial motion trajectory guess. This initial guess provides a simple starting point for motion, assuming that the valve moves at a constant speed in the sequence.

[0184] Next, Gaussian pyramid downsampling is performed on the ultrasound image sequence and the initial motion trajectory guess in the time dimension. Specifically, this includes low-pass filtering and interval sampling operations on both along the time axis to reduce their time resolution, thereby constructing a multi-layer time resolution data from coarse to fine, including the downsampled ultrasound image sequence and the downsampled control point state sequence. The coarsest layer has the fewest frames and the largest time interval, corresponding to the lowest time resolution. For example, a setting with three pyramid layers: the coarsest layer (L1) downsamples the original number of frames to about 1 / 4, the middle layer (L2) downsamples to about 1 / 2, and the finest layer (L3) is the original resolution. The Gaussian filter on the time axis can use a small kernel with a standard deviation of 1 to moderately smooth the time signal before reducing the sampling rate.

[0185] Then, an iterative optimization process is executed, starting from the coarsest time resolution layer and progressing layer by layer to the finest layer, specifically:

[0186] R1. At the current time resolution layer, based on the downsampled ultrasound image sequence corresponding to this layer and the control point state sequence of the current layer, construct and minimize the constructed global time domain energy function, and solve for the motion trajectory of the control point in the current layer. At each layer, use a numerical optimization algorithm (such as the L-BFGS algorithm) to minimize the energy function. Since the coarse resolution layer has fewer variables (time points), the optimization problem has a lower dimension, the computational cost is lower, and it is easier to avoid local minima caused by noise, thus quickly finding a coarse trajectory close to the global optimum.

[0187] R2. If the current layer is not the finest original temporal resolution layer, the control point motion trajectory optimized in step R1 is upsampled to the next layer with a finer temporal resolution through interpolation on the time axis, and this upsampled trajectory is used as the initial control point state sequence for the next layer optimization process. Cubic spline interpolation can be used to ensure the smoothness of the trajectory. The upsampling process transmits the "skeleton" information of the coarse layer optimization to the fine layer, providing an excellent starting point for the fine layer optimization and greatly reducing the search space.

[0188] R3. Set the next layer as the new current layer and repeat R1 and R2.

[0189] The optimization process starts from the coarsest temporal resolution layer and iterates through steps R1 to R3 until the global temporal energy function is minimized at the finest original temporal resolution layer. Finally, the physical elastic curve model is output as the control point motion trajectory over the entire time range of the original ultrasound image sequence. This "coarse-to-fine" strategy effectively solves the problem of high-dimensional nonlinear optimization easily getting trapped in local optima, and takes into account both the global convergence efficiency of the algorithm and the spatiotemporal accuracy of the final motion trajectory.

[0190] In this embodiment, it is particularly important to explain the process in the analysis module of decomposing the effective motion displacement from the control point motion trajectory and calculating the curvature-weighted time delay between the leaflets. The specific process is as follows:

[0191] First, based on the received control point motion trajectory, the motion trajectory of a cardiac stable reference point is determined. The cardiac stable reference point is obtained by automatically tracking a preset, relatively fixed anatomical landmark in the ultrasound image sequence, such as the center point of the atrioventricular ring or a stable attachment point of the ventricular wall. Next, the motion trajectory of the cardiac stable reference point is subjected to time low-pass filtering to extract its overall motion component that changes slowly over time. The cutoff frequency of the filter should be significantly lower than the fetal heart rate to separate the low-frequency overall motion of the heart, such as translation and rotation, while preserving the high-frequency opening and closing motion of the valves themselves.

[0192] Based on this overall motion component, for each moment and each control point in the physical elastic curve model, perform a motion decomposition operation:

[0193] S1. Calculate the absolute position change of the control point from the initial time to the current time. The change is the position coordinate of the control point at the current time minus its position coordinate at the initial time.

[0194] S2. Calculate the rigid motion displacement that the control point should undergo, induced by the overall motion of the cardiac stability reference point. This rigid motion displacement is obtained through the following steps:

[0195] S2.1 Obtain the relative position vector of the control point relative to the initial global motion component of the cardiac stability reference point at the initial moment. This vector is the position coordinate of the control point at the initial moment minus the coordinate of the global motion component of the reference point at the initial moment.

[0196] S2.2 Obtain an optimal rotation matrix, which is obtained by solving a minimization problem. The goal of this minimization problem is to minimize the sum of squares of the total difference between the current position of all control points relative to the reference point in the current overall motion component and the result obtained by applying the optimal rotation matrix to the relative position vector of each control point relative to the initial overall motion component of the reference point at the initial time. This step is equivalent to finding an optimal rotation at the current time that makes the configuration of all control points relative to the reference point as aligned as possible with the configuration at the initial time, thereby solving for the rotation component of the heart as a whole.

[0197] S2.3 Apply the optimal rotation matrix to the relative position vector obtained in step S2.1;

[0198] S2.4 Add the displacement vector of the overall motion component of the reference point from the initial time to the current time to the result of step S2.3. The displacement vector is the overall motion component coordinate of the reference point at the current time minus its overall motion component coordinate at the initial time.

[0199] S2.5. Subtract the position coordinates of the control point at the initial moment from the result obtained in step S2.4 to finally obtain the rigid motion displacement;

[0200] S3. Subtract the rigid motion displacement calculated in step S2 from the absolute position change calculated in step S1. The result is the effective motion displacement of the control point at that moment. Through this operation, the effective motion displacement filters out the "accompanying motion" caused by the movement and rotation of the heart as a whole in the chest cavity, thus reflecting the active deformation of the valve leaflets such as opening, closing, extension and bending more purely.

[0201] Through the above processes S1 to S3, the effective motion displacement of all control points at all times is obtained, forming an effective displacement vector field.

[0202] After obtaining the effective displacement vector field, curvature-weighted time delay is calculated. The process is as follows: For each control point, the rate of change of the magnitude of its effective motion displacement vector at adjacent time points is calculated, which is used as the effective motion rate of that control point. A rate curve is generated that changes over time. The time point corresponding to the peak of this rate curve within a complete cardiac cycle is identified and defined as the effective peak motion time of that control point. Simultaneously, based on the geometry of the physical elastic curve model at the initial time, the local curvature value at each control point is calculated. This local curvature value reflects the degree of curvature of the curve at that point; typically, the curvature value is highest in the valve tip region. For a pair of control points for which motion synchronicity needs to be evaluated, the effective peak motion time and its local curvature value for both control points are obtained. The local curvature value of each control point is multiplied by a curvature sensitivity coefficient greater than zero. Then, this product is used as the exponent of the natural constant e for calculation, and the result is the curvature weight of the control point. This calculation gives a significantly larger weight to high curvature points. This exponential mapping method can significantly amplify the contribution of the time information of high curvature regions (such as valve tips) to the time delay calculation, because the motion delay of the valve tip is usually more clinically significant than the delay of the valve body. The curvature sensitivity coefficient can be between 0.5 and 2.0, for example, it can be 1.0. Finally, the absolute value of the difference between the curvature weights of the two control points and their respective peak times is calculated, and this absolute value is divided by the sum of the maximum value of the two curvature weights and a preset normal coefficient. The result is defined as the curvature-weighted time delay between the two control points. This time delay value is used to quantify the degree of motion desynchronization between the two. The normal coefficient is used to prevent the denominator from being zero. The normal coefficient can be 0.001.

[0203] The motion vector field on the valve surface is obtained by interpolating the effective motion displacements of all control points, and the local strain rate is calculated. The specific process is as follows:

[0204] First, a continuous motion vector field is constructed based on the effective displacement vector field. The process is as follows: For each moment in the ultrasound image sequence, the static position coordinates of all control points of the physical elastic curve model at the initial moment are used as the source point for interpolation, and the effective motion displacement vectors corresponding to each control point at the same moment are used as the target value of the source point. A two-dimensional vector field is fitted on all source points using the thin-plate spline interpolation method. Specifically, this is achieved by minimizing a target function consisting of a data fitting term and a regularization term. The data fitting term is calculated as the sum of the squares of the Euclidean norms of the differences between the two-dimensional vector field function values ​​to be fitted at all source points and the target effective displacement vectors corresponding to those points. The regularization term is used to penalize excessive curvature of the vector field function, and its value is a preset regularization parameter multiplied by an integral expression over the entire two-dimensional image plane. The calculation result, specifically the integral expression, is as follows: the square of the second-order partial derivative of the component function of the vector field function to be fitted with respect to the x-coordinate, plus twice the square of the mixed second-order partial derivative of the component function of the x-coordinate with respect to the x and y coordinates, plus the square of the second-order partial derivative of the component function of the x-coordinate with respect to the y-coordinate, plus twice the square of the mixed second-order partial derivative of the component function of the y-coordinate with respect to the x and y coordinates, plus the square of the second-order partial derivative of the component function of the y-coordinate with respect to the y-coordinate, and finally summing these six results and integrating them over the entire plane; this regularization term is essentially a measure of the total bending energy of the vector field, and its function is to ensure that the vector field interpolated between control points changes smoothly, avoiding drastic physical fluctuations. The regularization parameter controls the smoothness, and its value needs to be determined during debugging based on specific data. By solving the problem of minimizing the objective function, a continuous motion vector function defined on the two-dimensional image plane is obtained. The value of this function at any point represents the effective motion velocity vector of that point at the current moment, thus obtaining the continuous motion vector field of the valve surface.

[0205] After obtaining the continuous motion vector field, the local strain rate is calculated. The process is as follows: For each location point within the continuous motion vector field, the following steps are performed sequentially:

[0206] T1. Calculate the spatial velocity gradient tensor at this location: Decompose the velocity of the continuous motion vector field at the current point into horizontal and vertical components; calculate the partial derivatives of the horizontal velocity component with respect to the horizontal coordinate and the vertical coordinate, and the partial derivatives of the vertical velocity component with respect to the horizontal coordinate and the vertical coordinate, respectively; these four partial derivative values ​​together form a 2x2 matrix, called the spatial velocity gradient tensor at this point; the partial derivatives are approximated on a regular grid using the central difference method;

[0207] T2. Calculate the strain rate tensor at this location: Add the spatial velocity gradient tensor obtained in step T1 to its own transpose matrix, and then multiply the sum by one-half. The result is called the strain rate tensor at this point. This tensor is a symmetric second-order tensor, and its components represent the local tensile or compressive rates of the medium at this point in different directions. The strain rate tensor is the core physical quantity describing the local deformation rate in continuum mechanics.

[0208] T3. Extract one or more scalarized local strain rate indices from the strain rate tensor, specifically including:

[0209] T3.1 Maximum Principal Strain Rate: The eigenvalues ​​of the strain rate tensor obtained in step T2 are used as the maximum principal strain rate at that point. This value characterizes the rate of maximum elongation or shortening in the local deformation at that point. This index is particularly sensitive for identifying areas on the valve that are excessively stretched (which may lead to weakness) or abnormally compressed.

[0210] T3.2 Equivalent tensile-compressive strain rate: The sum of the squares of the components of the strain rate tensor obtained in step T2 is multiplied by two-thirds, and then the square root of the product is taken as the equivalent tensile-compressive strain rate at that point. This value comprehensively reflects the overall deformation rate intensity caused by pure tension or compression at that point. "Sum of the squares of the components" specifically refers to the sum of the squares of the values ​​of all four components of the strain rate tensor (usually three independent components due to symmetry).

[0211] T3.3 Shear strain rate: Calculate the determinant of the strain rate tensor obtained in step T2 and calculate the square of its trace; when the trace of the strain rate tensor is not zero, divide four times the determinant value by the square of the trace, and then take the square root of the quotient. The result is taken as the shear strain rate at that point. This value characterizes the rate of local shape change caused by shear deformation at that point; the shear strain rate helps to identify sliding or misalignment motion modes between or within the leaflets.

[0212] By sequentially performing steps T1 to T3 on all locations on the valve surface, the maximum principal strain rate distribution map, the equivalent tensile-compressive strain rate distribution map, and the shear strain rate distribution map covering the entire valve region are obtained. These distribution maps together constitute the local strain rate, providing clinicians with quantitative and visual information on the local mechanical load state of the valve tissue, and assisting in the judgment of its functional abnormalities and potential risks.

[0213] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0214] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0215] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0216] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0217] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0218] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0219] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A system for identifying and assessing valvular motion abnormalities in fetal echocardiography, characterized in that, Specifically, it includes: It includes an initialization module, a modeling module, an optimization module, and an analysis module that are sequentially connected, among which; Initialization module: Used to receive ultrasound image sequence of fetal heart, select a frame image from the ultrasound image sequence, perform morphological opening operation on a set of linear structuring elements with increasing scale for the single frame image, calculate the morphological gradient of the image after opening operation at each scale, preserve the ridge structure that is stable at multiple consecutive scales through scale space analysis, and generate and output a skeleton with a single pixel width by performing a final erosion process that preserves the connection of the union of the preserved stable ridge structures. Modeling module: It receives a skeleton with a width of one pixel, detects the branch points and endpoints of the skeleton to use them as fixed key control points, and adaptively inserts additional control points on the skeleton segments between the key points according to the local curvature changes, thereby generating a set of control points. Based on this set of control points, it constructs a physical elastic curve model containing spring damper connections and discrete bending energy terms to complete the initialization of the physical elastic curve model. The optimization module receives the ultrasound image sequence and the initialized physical elastic curve model. It constructs a global time-domain energy function and minimizes the function using a multi-resolution optimization strategy to obtain the control point motion trajectory of the physical elastic curve model over the entire image sequence time range. The global time-domain energy function includes data terms, internal constraint terms, and time-domain smoothing terms. Analysis module: It is used to receive the motion trajectory of the control point, decompose the effective motion displacement relative to the cardiac stability reference point from the motion trajectory of the control point, and calculate the curvature weighted time delay between the leaflets based on the time when the effective motion of each control point reaches the peak and the local curvature of its motion trajectory to quantify motion desynchronization. At the same time, it obtains the motion vector field of the valve surface by interpolation based on the effective motion displacement of all control points and calculates the local strain rate.

2. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 1, characterized in that: In the initialization module, a set of linear structuring elements with increasing scale are used to perform morphological opening operations sequentially, and the morphological gradient of the image after the opening operation at each scale is calculated. Specifically, this includes: Select a frame from the received ultrasound image sequence and define a set of scale-increasing linear structuring elements; For each structural element at each scale, a morphological opening operation is first performed on the selected image to obtain the opening operation image at that scale. Subsequently, for each opening operation image, a first intermediate image is obtained by performing a dilation operation using a preset compact structuring element with a size smaller than all linear structuring elements, and a second intermediate image is obtained by performing an erosion operation. The values ​​of corresponding pixels in the first intermediate image and the second intermediate image are subtracted to obtain the morphological gradient image corresponding to that scale, thereby generating a set of morphological gradient images that correspond one-to-one with multiple scales. Preserving ridge structures that are stable across multiple consecutive scales through scale-space analysis involves the following steps: For a set of morphological gradient images, for each pixel position in a frame of an image, the following judgment and calculation are performed: determine the scale at which the gradient value of the pixel position first exceeds a preset gradient threshold in a set of morphological gradient images, and take it as the starting scale; starting from the starting scale, along the scale increment direction, count the total number of adjacent scales in which the gradient value of the pixel position continuously exceeds the gradient threshold, and define the total number as the scale persistence measure of the pixel position. Simultaneously, a set of gradient values ​​is obtained at these continuous scales where the gradient value at the pixel location continuously exceeds the gradient threshold. The maximum and minimum values ​​in this set of gradient values ​​are found, and the difference between the maximum and minimum values ​​is calculated to obtain the range of the set of gradient values. Then, the ratio of the range to a corrected denominator is calculated. The corrected denominator is the sum of the maximum value in the set of gradient values ​​and a preset minimum positive constant. Finally, the ratio is subtracted from the number 1, and the result is defined as the gradient consistency measure at the pixel location. A minimum continuous scale number threshold and a gradient consistency threshold are preset. For each pixel position in a frame of image, the pixel is determined to belong to a ridge structure that is stable at multiple consecutive scales only if the scale persistence metric of the pixel position is not less than the minimum continuous scale number threshold and the gradient consistency metric of the pixel position is higher than the gradient consistency threshold. All pixel locations identified as belonging to stable ridge structures are marked to generate a binary stable ridge marker image.

3. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 2, characterized in that: After generating the stable ridge marker image, the union of the preserved stable ridge structures is subjected to a final erosion process that preserves connectivity. The specific process is as follows: First, morphological closing operations are performed on the stable ridge marker image to connect the small gaps in the image caused by noise or insufficient contrast, resulting in a aggregated ridge image. Calculate the morphological skeleton of the aggregated ridge image to obtain a preliminary skeleton image composed of single-pixel wide foreground pixels; Calculate the distance transformation map of the aggregated ridge image. The value of each foreground pixel in the distance transformation map represents the distance from that point to the nearest background pixel in the image. Detect skeleton branch point pixels in the preliminary skeleton image. Skeleton branch point pixels refer to pixels in the preliminary skeleton image that are adjacent to three or more other skeleton pixels. Based on the distance transformation map and the positions of the detected skeleton branch point pixels, a connectivity-preserving iterative morphological erosion operation is performed on the aggregated ridge image. The rules of the iterative morphological erosion operation are as follows: In each iteration, only foreground pixels in the aggregated ridge image that simultaneously meet the following two conditions are allowed to be removed: First, the pixel is not located within a protected neighborhood centered on any skeleton branch point pixel with a preset pixel distance as its radius; Second, the distance value corresponding to the pixel in the distance transformation map is greater than the number of erosion iterations already performed. Repeatedly perform the iterative morphological erosion operation until the width of all remaining foreground connected regions in the aggregated ridge image becomes one pixel, or until a preset maximum number of iterations is reached. At this point, stop the iteration and use the resulting image as a single-pixel-wide skeleton.

4. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 3, characterized in that: After obtaining the single-pixel-wide skeleton, a skeleton post-processing operation is included before output, which specifically includes: Branch analysis and cleanup are performed on the single-pixel wide skeleton image obtained after the final erosion process that preserves connectivity. The branch analysis and cleanup are achieved through the following steps: identifying all endpoint pixels in the single-pixel wide skeleton image, where an endpoint pixel is a pixel in the single-pixel wide skeleton image that is only adjacent to one other skeleton pixel. For each identified endpoint pixel, starting from that endpoint pixel, traverse along the path of the skeleton pixels connected to it, and record all the pixels traversed until a skeleton branch point pixel or another endpoint pixel is encountered. The skeleton branch point pixel refers to a pixel in the single-pixel wide skeleton image that is connected to three or more other skeleton pixels. The total number of pixels traversed from the starting endpoint pixel to the ending point is defined as the length of the skeleton branch originating from that starting endpoint pixel. A preset length threshold is set; if the calculated length of a certain skeleton branch is less than the length threshold, the skeleton branch is determined to be a short branch caused by noise, and all pixels constituting the skeleton branch are deleted from the single-pixel wide skeleton image. After performing the above traversal, length calculation, judgment and pixel deletion operations on all identified endpoint pixels, the cleaned skeleton containing only the main trunk and significant branches with a width of one pixel is finally obtained and output.

5. The fetal cardiac ultrasound system for identifying and assessing valvular motion abnormalities in fetuses according to claim 4, characterized in that: In the modeling module, branch points and endpoints of the skeleton are detected and used as fixed key control points. Additional control points are then adaptively inserted into the skeleton segments between these key points based on local curvature changes. The specific process is as follows: First, traverse the received single-pixel-wide skeleton, identify all skeleton branch point pixels and endpoint pixels, and define all identified skeleton branch point pixels and endpoint pixels together as a set of fixed key control points. Subsequently, on the skeleton path segment determined by any two adjacent fixed key control points, additional control points are adaptively inserted based on the local curvature values ​​of each pixel on the path segment. The calculation process of the local curvature value is as follows: For a pixel to be calculated on the path segment, obtain the position coordinates of its direct predecessor pixel and direct successor pixel on the path segment, calculate the first vector from the predecessor pixel to the pixel to be calculated, and the second vector from the pixel to be calculated to the successor pixel, and calculate the angle between the first vector and the second vector; then calculate the sum of the magnitudes of the first vector and the second vector; finally, multiply the value by the absolute value of the sine of the angle, and divide by the sum to obtain the local curvature value of the pixel to be calculated. After calculating the local curvature values ​​of each pixel on the skeleton path segment, the local curvature values ​​are compared with the preset high curvature threshold and medium curvature threshold. The comparison process is as follows: in the path segment where the local curvature value is greater than the high curvature threshold, additional control points are inserted at the first preset interval distance. In path segments where the local curvature value is between the medium curvature threshold and the high curvature threshold, additional control points are inserted at a second preset interval distance that is greater than the first preset interval distance. In path segments where the local curvature value is less than the medium curvature threshold, insert additional control points at a third preset interval distance that is greater than the second preset interval distance, or do not insert additional control points. The additional control points generated through the above comparison process constitute an additional control point set. The union of the fixed key control point set and the additional control point set is the generated set of control points. The initial position of each control point is set to the coordinates of the corresponding pixel in the single-pixel-wide skeleton.

6. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 5, characterized in that: Based on this set of control points, a physical elastic curve model including spring-damper connections and discrete bending energy terms is constructed. The specific process is as follows: Each point in a set of control points is defined as a mass point, and adjacent points on a single-pixel-wide skeleton in this set of control points are connected using spring-damper units to form a chain structure. After the connection is completed, based on the geometric information of the single-pixel-wide skeleton, the parameters in the physical elastic curve model are initialized. The initialization includes: Q1. For each pair of adjacent control points connected by spring damper units, define a natural length parameter, the value of which is set to the shortest path pixel distance between the two control points in a single-pixel-wide skeleton multiplied by a preset pixel physical size. Q2. For each pair of adjacent control points connected by spring damper units, define a spring stiffness coefficient, which is set to be proportional to the reciprocal of the local or average width of the skeleton in the corresponding region. Q3. Define a bending stiffness coefficient for any three consecutive control points in a chain structure. For control points located near the skeleton branch point pixels in the fixed key control point set, the associated bending stiffness coefficient is set to a higher value than that of the non-branch point region. Q4. For any three control points arranged consecutively in a chain structure, define a balance angle parameter. Its value is directly set to the exterior angle value calculated by the three control points based on their initial position coordinates in a single-pixel-wide skeleton. Q5. For each control point or each spring damper connection, define a damping coefficient that is proportional to the product of the square root of the corresponding spring stiffness coefficient and the mass parameter assigned to the control point. After completing the above initialization assignment, define the internal elastic potential energy of the physical elastic curve model. The internal elastic potential energy consists of two parts: spring tension energy and discrete bending energy. By defining the internal elastic potential energy and initializing the parameters as described above, the physical elastic curve model is initialized.

7. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 6, characterized in that: The optimization module constructs a global time-domain energy function, and the specific process is as follows: First, the optimization module receives the ultrasound image sequence and the initialized physical elastic curve model. The state of the physical elastic curve model at any time is defined by the position of all its control points at that time. The solution objective is to obtain the motion trajectory of all control points of the physical elastic curve model over the entire time range of the ultrasound image sequence. To solve for the motion trajectory, a global time-domain energy function is constructed. This function is a weighted sum of the data terms, internal constraint terms, and time-domain smoothing terms. Specifically, the global time-domain energy function equals the sum of the following three terms: The first step is to calculate the data items for each moment in the ultrasound image sequence from the first moment to the last moment, and then sum all the data items. The second term involves multiplying an internal constraint term weight coefficient by the internal constraint term for each moment in the ultrasound image sequence from the first moment to the last moment, calculating the internal constraint term for that moment, and summing all these internal constraint terms. The third term is a time-domain smoothing term weighting coefficient multiplied by each time point from the second time point to the last time point in the ultrasound image sequence, calculating the time-domain smoothing term between that time point and the previous time point, and summing all these time-domain smoothing terms. The calculation process for the data items is as follows: First, for each frame of the ultrasound image sequence, calculate its phase consistency feature map and local principal direction field; Then, for the continuous curve defined by all control point positions at a certain moment in the physical elasticity curve model, a series of points are sampled on the curve; for each sampled point, the following operations are performed: The feature intensity value at the sampling point is extracted from the phase consistency feature map. Based on the local shape of the continuous curve near the sampling point, the tangent direction of the curve at the point is calculated and normalized into a vector of unit length, called the unit tangent vector of the sampling point. The local principal direction angle value of the image structure at the sampling point is extracted from the local principal direction field and converted into a unit direction vector. The absolute value of the angle between the unit tangent vector and the unit direction vector of the sampling point is calculated. A Gaussian kernel function with a preset standard deviation parameter is applied to the absolute value to map it, and a direction consistency weight value between zero and one is obtained. Finally, the feature intensity value is multiplied by the orientation consistency weight value to obtain the orientation matching contribution value of the sampling point; the orientation matching contribution values ​​of all sampling points are summed and negatively taken as the main part of the data item at that moment. This data item also includes a distance penalty sub-item. The final value of the data item is obtained by adding the sum of the negative orientation matching contribution values ​​to the value of the distance penalty sub-item. The calculation process of the internal constraint term is as follows: at each moment, the internal elastic potential energy of the physical elastic curve model at that moment is directly calculated. The internal elastic potential energy is composed of the spring tension energy and the discrete bending energy. The value of the internal constraint term is equal to the sum of the spring tension energy and the discrete bending energy. The calculation process of the time-domain smoothing term is as follows: For any two adjacent time points, denoted as the previous time point and the current time point, the time-domain smoothing term is composed of a weighted sum of two parts: a smoothness penalty term based on motion acceleration and a penalty term based on motion pattern consistency. The calculation process for the smoothness penalty term based on motion acceleration is as follows: For each control point in the physical elastic curve model, firstly, calculate its displacement vector from the previous moment to the current moment. This displacement vector is equal to the control point's position coordinate at the current moment minus its position coordinate at the previous moment. Simultaneously, calculate the control point's displacement vector relative to the even earlier moment from the previous moment. Then, calculate the difference between the two displacement vectors to obtain an approximate acceleration vector of the control point at the current moment. Next, calculate the square of the magnitude of this approximate acceleration vector, multiply it by a first acceleration weighting coefficient, and obtain the first acceleration weighting coefficient of the control point. First, an acceleration penalty value is calculated. Simultaneously, the ratio of the approximate magnitude of the acceleration vector to a correction denominator is calculated, and this ratio is squared and multiplied by a second acceleration weighting coefficient to obtain the second acceleration penalty value for that control point. The correction denominator is the sum of the magnitude of the control point's displacement vector at the previous moment and 0.

001. Finally, a preset control point weighting coefficient, the first acceleration penalty value, and the second acceleration penalty value are multiplied together to obtain the acceleration smoothing penalty contribution value for that control point. The acceleration smoothing penalty contribution values ​​for all control points are summed to obtain the smoothness penalty term based on motion acceleration. The calculation process for the motion pattern consistency-based penalty term is as follows: First, the sets of displacement vectors of all control points at the previous time step and the sets of displacement vectors of all control points at the current time step are statistically analyzed. Each set is considered as an empirical probability distribution of a displacement vector. Then, the Kullback-Leibler divergence of the empirical distribution of the displacement vector at the current time step relative to the empirical distribution of the displacement vector at the previous time step is calculated. Finally, the Kullback-Leibler divergence is multiplied by a pattern consistency weight coefficient to obtain the motion pattern consistency-based penalty term. The value of the time-domain smoothing term is equal to the sum of the smoothness penalty term based on motion acceleration and the penalty term based on motion pattern consistency.

8. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 7, characterized in that: The multi-resolution optimization strategy specifically includes: First, based on the initial static positions of each control point in the initialized physical elastic curve model received from the modeling module, a preliminary control point position sequence that is uniformly distributed in time and runs through the entire ultrasound image sequence is generated using a linear interpolation method, and this preliminary sequence is defined as the initial motion trajectory guess. Next, Gaussian pyramid downsampling is performed on the ultrasound image sequence and the initial motion trajectory guess in the time dimension. Specifically, low-pass filtering and interval sampling are performed on both along the time axis to construct a multi-layer time resolution data from coarse to fine, including the downsampled ultrasound image sequence and the downsampled control point state sequence. Then, an iterative optimization process is executed, starting from the coarsest time resolution layer and progressing layer by layer to the finest layer, specifically: R1. At the current time resolution layer, based on the downsampled ultrasound image sequence corresponding to this layer and the control point state sequence of the current layer, construct and minimize the constructed global time domain energy function, and solve to obtain the motion trajectory of the control point of the current layer. R2. If the current layer is not the finest original temporal resolution layer, the control point motion trajectory obtained in step R1 is upsampled to the next layer with a finer temporal resolution through interpolation on the time axis, and this upsampled trajectory is used as the initial control point state sequence of the next layer optimization process. R3. Set the next layer as the new current layer and repeat R1 and R2. The optimization process starts from the coarsest temporal resolution layer and iterates through steps R1 to R3 until the global temporal energy function is minimized at the finest original temporal resolution layer. Finally, the physical elastic curve model is output as the control point motion trajectory over the entire time range of the original ultrasound image sequence.

9. The fetal cardiac ultrasound system for identifying and assessing valvular motion abnormalities in fetuses according to claim 8, characterized in that: In the analysis module, the effective motion displacement is decomposed from the motion trajectory of the control point and the curvature-weighted time delay between the leaflets is calculated. The specific process is as follows: First, based on the received motion trajectory of the control point, the motion trajectory of a cardiac stability reference point is determined; the cardiac stability reference point is obtained by automatically tracking a preset, relatively fixed anatomical landmark in the ultrasound image sequence; then, the motion trajectory of the cardiac stability reference point is subjected to time low-pass filtering to extract its overall motion component that changes slowly over time. Based on this overall motion component, for each moment and each control point in the physical elastic curve model, a motion decomposition operation is performed: S1. Calculate the absolute position change of the control point from the initial time to the current time. The change is the position coordinate of the control point at the current time minus its position coordinate at the initial time. S2. Calculate a rigid motion displacement that the control point should undergo, induced by the overall motion of the cardiac stability reference point. This rigid motion displacement is obtained through the following steps: S2.1 Obtain the relative position vector of the control point relative to the initial overall motion component of the cardiac stability reference point at the initial moment. The vector is the position coordinate of the control point at the initial moment minus the coordinate of the overall motion component of the reference point at the initial moment. S2.2 Obtain an optimal rotation matrix, which is obtained by solving a minimization problem. The goal of the minimization problem is to minimize the sum of squares of the total difference between the position of all control points relative to the current overall motion component of the reference point at the current time and the result obtained by applying the optimal rotation matrix to the relative position vector of each control point relative to the initial overall motion component of the reference point at the initial time. S2.3 Apply the optimal rotation matrix to the relative position vector obtained in step S2.1; S2.4 Add the result of step S2.3 to the displacement vector of the overall motion component of the reference point from the initial time to the current time. The displacement vector is the overall motion component coordinate of the reference point at the current time minus its overall motion component coordinate at the initial time. S2.

5. Subtract the position coordinates of the control point at the initial moment from the result obtained in step S2.4 to finally obtain the rigid motion displacement; S3. Subtract the rigid motion displacement calculated in step S2 from the absolute position change calculated in step S1. The result is the effective motion displacement of the control point at that moment. Through the above processes S1 to S3, the effective motion displacement of all control points at all times is obtained, forming an effective displacement vector field. After obtaining the effective displacement vector field, the curvature-weighted time delay is calculated.

10. The fetal cardiac ultrasound system for identifying valvular motion abnormalities and assessing function according to claim 9, characterized in that: The specific process of obtaining the motion vector field of the valve surface and calculating the local strain rate based on the effective motion displacement interpolation of all control points is as follows: First, a continuous motion vector field is constructed based on the effective displacement vector field. The process is as follows: For each moment in the ultrasound image sequence, the static position coordinates of all control points of the physical elastic curve model at the initial moment are used as the source point for interpolation, and the effective motion displacement vectors corresponding to each control point at the same moment are used as the target value of the source point. A two-dimensional vector field is fitted on all source points using the thin-plate spline interpolation method. Specifically, this is achieved by minimizing an objective function consisting of a data fitting term and a regularization term. The data fitting term is calculated as the sum of the squares of the Euclidean norms of the differences between the two-dimensional vector field function values ​​to be fitted at all source points and the target effective displacement vectors corresponding to those points. The regularization term is used to penalize excessive curvature of the vector field function. Its value is a preset regularization parameter multiplied by the result of an integral expression calculated over the entire two-dimensional image plane. By solving the problem of minimizing the objective function, a continuous motion vector function defined on the two-dimensional image plane is obtained. The value of this function at any point represents the effective motion velocity vector of that point at the current moment, thus obtaining the continuous motion vector field of the valve surface. After obtaining the continuous motion vector field, the local strain rate is calculated as follows: For each location point within the continuous motion vector field, the following steps are performed sequentially: T1. Calculate the spatial velocity gradient tensor at this location: Decompose the velocity of the continuous motion vector field at the current point into horizontal and vertical components; calculate the partial derivatives of the horizontal velocity component with respect to the horizontal coordinate and the vertical coordinate, respectively; these four partial derivative values ​​together form a 2x2 matrix, called the spatial velocity gradient tensor at this point; the partial derivatives are approximated on a regular grid using the central difference method; T2. Calculate the strain rate tensor at this location: Add the spatial velocity gradient tensor obtained in step T1 to its own transpose matrix, and then multiply the result by one-half. The result is called the strain rate tensor at this point. This tensor is a symmetric second-order tensor, and its components represent the local tensile or compressive rates of the medium at this point in different directions. T3. Extract one or more scalarized local strain rate indices from the strain rate tensor, specifically including: T3.1 Maximum Principal Strain Rate: The eigenvalues ​​of the strain rate tensor obtained in step T2 are used as the maximum principal strain rate at that point. This value characterizes the rate of maximum elongation or shortening in the local deformation at that point. T3.2 Equivalent tensile and compressive strain rate: The sum of the squares of the components of the strain rate tensor obtained in step T2 is multiplied by two-thirds, and then the square root of the product is taken as the equivalent tensile and compressive strain rate at that point. This value comprehensively reflects the overall deformation rate intensity caused by pure tension or compression at that point. T3.3 Shear strain rate: Calculate the determinant value of the strain rate tensor obtained in step T2, and calculate the square of its trace; when the trace of the strain rate tensor is not zero, divide four times the determinant value by the square of the trace, and then take the square root of the quotient value. The result is taken as the shear strain rate at that point. By sequentially performing steps T1 to T3 on all locations on the valve surface, the maximum principal strain rate distribution map, the equivalent tensile-compressive strain rate distribution map, and the shear strain rate distribution map covering the entire valve region are obtained. These distribution maps together constitute the local strain rate.

Images