Heart motion feature extraction method based on optical flow estimation

Abstract

The invention discloses a heart motion feature extraction method based on optical flow estimation, and relates to medical image processing. Preprocessing the input four-dimensional space-time cardiac magnetic resonance imaging data, namely, scaling pixel values; inputting two frames of images which are continuous in time into a feature encoder and a context encoder for feature extraction, wherein the two frames of images are divided into a reference frame and a moving frame; calculating the correlation between the feature maps through a correlation volume calculation module, constructing a correlation pyramid, and extracting the feature maps to provide matching information for subsequent optical flow estimation; the motion feature iteration enhancement module iteratively and continuously refines an optical flow estimation result through a deformable convolution and global motion aggregation (GMA) module; model parameters are optimized based on a weighted sum of luminosity consistency loss, smoothness loss, and gradient consistency loss. Edge features are adaptively captured through deformable convolution, global and local features are fused by using a GMA module, optical flow prediction errors are effectively reduced, motion estimation quality is improved, and key features of a heart edge region are maintained.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and in particular to an improved method for extracting cardiac motion features based on optical flow estimation by Recurrent All-Pairs Field Transforms with Global Motion Aggregation (RAFT) and Global Motion Aggression (GMA), which is applicable to the motion analysis of cardiac magnetic resonance imaging (MRI) data. Background Art

[0002] The accurate extraction of cardiac motion features is crucial for the diagnosis and functional evaluation of heart diseases. Traditional cardiac motion analysis methods mainly rely on manual measurements based on specific anatomical landmark points or analysis methods based on deformation tracking. These methods have problems such as complex operations and being easily affected by image quality, and it is difficult to meet the efficient and accurate clinical needs. In recent years, deep learning-based optical flow estimation algorithms have been widely used in cardiac motion extraction. Optical flow describes the motion information of pixels between adjacent time frames and can effectively depict the dynamic changes of the myocardium. With the development of Convolutional Neural Networks (CNNs), many deep learning models have been proposed for optical flow estimation in cardiac magnetic resonance images. PWC-Net (Sun D, Yang X, Liu M-Y, et al. PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume [C]. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018: 8934-8943) combines a pyramid model and a cost volume to refine the optical flow estimation layer by layer in a multi-scale space and can efficiently process large displacement scenarios. RAFT (Teed Z, Deng J. RAFT: recurrent all-pairs field transforms for optical flow [C], 2020: 402-419) iteratively updates and refines the optical flow by extracting features and constructing a four-dimensional volume to obtain a high-precision optical flow field.

[0003] Currently, there are still the following two problems in the deep learning-based cardiac motion feature extraction model: (1) In the case of large cardiac motion, the accuracy of motion estimation is insufficient; (2) Due to the fact that the optical flow method usually adopts a smoothness constraint, in the edge region of the cardiac structure, the motion estimation error is large. Summary of the Invention

[0004] The object of the present invention is to provide a method for extracting cardiac motion features based on optical flow estimation to solve the problem of insufficient estimation accuracy of existing methods in large motion areas and edge areas, and to achieve more accurate extraction of cardiac motion features through a deep learning framework.

[0005] To achieve the above object of the invention, the present invention provides the following technical solutions.

[0006] The present invention provides a method for extracting cardiac motion features based on optical flow estimation, including the following steps:

[0007] 1) Preprocess the input four-dimensional spatio-temporal cardiac magnetic resonance imaging data, and scale the pixel values;

[0008] 2) Divide two temporally consecutive frames of images into a reference frame and a moving frame, and input them into a feature encoding module for feature extraction. The feature encoding module includes a feature encoder and a context encoder, where the feature encoder encodes both frames of images, and the context encoder only encodes the reference frame;

[0009] 3) Calculate the correlation between feature maps through a correlation volume calculation module, construct a correlation pyramid and extract feature maps from it to provide matching information for subsequent optical flow estimation;

[0010] 4) The motion feature iterative enhancement module iteratively refines the optical flow estimation result through deformable convolution and a global motion aggregation (GMA) module;

[0011] 5) Optimize the model parameters based on the weighted sum of photometric consistency loss, smoothness loss, and gradient consistency loss.

[0012] In step 1), the preprocessing of the input four-dimensional spatio-temporal cardiac magnetic resonance imaging data, that is, normalization processing, and scaling the pixel values. The process of normalization processing can be expressed as:

[0013]

[0014] where X is the input four-dimensional spatio-temporal cardiac magnetic resonance data, with a size of A×B×S×fr, A is the height of the image, B is the width of the image, S is the number of slices of the image, and fr is the number of time frames of the image; X norm is the data obtained after normalization processing, X min and X max are the minimum and maximum values of the pixels in the same image; crop the two-dimensional image with a size of A×B to obtain an image with a height and width of H×W.

[0015] In step 2), the specific steps for feature extraction may be as follows:

[0016] The feature encoder adopts a convolutional layer structure to encode the two input frames of images into low-dimensional feature maps respectively; the encoder g θ (·) encodes the reference frame I1 and the moving frame I2 into features that are respectively scaled down by times in length and width, that is λ where θ is a learnable network parameter, H and W are respectively the height and width of the feature map, ch is the number of channels, and a normalization layer and a rectified linear unit are connected after the convolutional layer; the feature encoder further includes J residual modules, each residual module contains two convolutional layers, and each convolutional layer is followed by a normalization layer and a rectified linear unit; the feature encoding module further includes a context encoder h θ (·), λ is a learnable network parameter, the architecture of the context encoder is similar to that of the feature encoder, extracts features from the reference frame I1, is independent of the encoder g θ (·) respectively, and the output feature extraction results are respectively input into the correlation volume calculation module and the attention module.

[0017] In step 3), calculating the correlation between feature maps through the correlation volume calculation module, constructing a correlation pyramid and extracting feature maps from it to provide matching information for subsequent optical flow estimation, the specific steps may be as follows:

[0018] Performing a dot product on each corresponding pixel point in all the extracted feature pairs g θ (I1) and g θ (I2) can obtain the correlation volume C, and this process can be expressed as:

[0019]

[0020] where g θ (I1) and g θ (I2) are the extracted feature pairs, H and W are respectively the height and width of the feature map, and i, j, k, l, and h are coefficients for traversing their respective dimensions;

[0021] The correlation volume calculation module consists of two parts: a correlation pyramid and a correlation search; in the correlation pyramid, four pooling layers are used to pool the last two dimensions of the correlation volume C to construct a four-layer correlation pyramid {C 1 , C 2 , C 3 , C 4}; where the size of the Kth layer (K = 1, 2, 3, 4) of the correlation pyramid C K is In the correlation search operation L c in, L cExtract feature maps from the correlation pyramid; given the current optical flow estimate (f 1 ,f 2 ), estimate the corresponding pixel U′=(u+f 1 (u),v+f 2 (v)); Define the search domain N(U′) of U′ r

[0022]

[0023] Among them, r is the search radius, dU is the offset, which belongs to the two-dimensional integer space, and ||·||1 is the L1 norm; for the relevant volume C k , in the search domain N(U′ / 2 K ) r Find the feature map (K=1,2,3,4) in the dataset, and finally concatenate the four feature maps into a correlation tensor.

[0024] In step 4), the motion feature iterative enhancement module iteratively refines the optical flow estimation result through deformable convolution and global motion aggregation GMA (Global Motion Aggression) modules, and the specific steps are as follows:

[0025] The motion feature iterative enhancement module inputs the initial estimate of optical flow and the context features processed by the attention module, and performs a deformable convolution operation on the two. In the deformable convolution operation, a regular grid is sampled on the input feature map. Increased offset {Δp m |m=1,...,M}, where p m Represents traversal The position in the convolution kernel, M is the number of sampling points in the convolution kernel, that is, |·| is the number of elements in the set; for each position p0 on the output feature map Y, there is:

[0026]

[0027] in, For a regular network, Δp m is the offset, p m Represents traversal of a regular grid The position in the m ) is the convolution kernel at p m The value at , T is the input feature map;

[0028] The input feature map T is implemented by bilinear interpolation as follows:

[0029]

[0030] where p represents any position of the input feature map (p = p0 + p m + Δp m ), q traverses all integer spatial positions in the input feature map T, and G(·, ·) is a bilinear interpolation kernel;

[0031] The current optical flow information obtained after the deformable convolution operation forms a joint representation of the motion information and the spatial matching relationship with the correlation volume, and global weighted modulation is applied to this joint representation through Global Motion Aggression (GMA); in the Global Motion Aggression (GMA) module, given the context feature and the motion feature where N = H × W and W are the height and width of the feature map respectively; represents the i-th context feature vector, represents the i-th motion feature vector, D c and D m are the numbers of the context feature vector and the motion feature vector respectively; the aggregated motion feature can be expressed as:

[0032]

[0033] where, is the updated i-th motion feature vector, represents the i-th motion feature vector, α is a learnable scalar parameter, φ and σ are the projection functions of the query vector, key vector, and value vector respectively; f is a similarity attention function, which can be expressed as:

[0034]

[0035] where D is the number of channels of the feature map, φ and σ are the projection functions of the query vector, key vector, and value vector respectively, and are expressed as:

[0036]

[0037] φ(x i ) = W key x i

[0038] σ(x i ) = W val x i

[0039] where, is the weight matrix, Denote the i-th context feature vector;

[0040] The output of the Global Motion Aggression (GMA) module is the concatenation of the context feature, the motion feature, and the aggregated motion feature After introducing the two-dimensional relative position embedding, the aggregated motion vector can be expressed as:

[0041]

[0042] Denote the i-th motion feature vector, α is a learnable scalar parameter, f is the similarity attention function, φ and σ are the projection functions of the query vector, key vector, and value vector respectively, z j-i denotes the relative position embedding vector of positions j and i; the model learns the vertical and horizontal offset embedding vectors respectively and adds them to obtain z j-i ;

[0043] In step 5), the model parameters are optimized based on the weighted sum of the photometric consistency loss, smoothness loss, and gradient consistency loss. The specific loss function is as follows:

[0044] For a pair of temporally consecutive image reference frame I1 and moving frame I2, the optical flow field v is used to reverse-deform I1 to obtain the image I1'; the photometric consistency loss L photo is defined as follows:

[0045]

[0046] where, is the value of the i-th pixel in the estimated image obtained by reverse deformation, is the value of the i-th pixel in the moving frame, N is the total number of image pixels, is the square of the L2 norm; the photometric consistency loss is used to measure the difference between the optical flow field deformed image and the real image. The smaller the photometric consistency loss, the more accurate the optical flow estimation;

[0047] The smoothness loss L smooth is defined as follows:

[0048]

[0049] where N is the total number of image pixels, v pred is the estimated optical flow field, is the optical flow vector of the i-th pixel of the estimated optical flow field, including the horizontal component and the vertical component; and It represents calculating the gradient of the optical flow field in the x-axis and y-axis directions; the smoothness loss is used to limit the discontinuity in the optical flow field, making the change of the optical flow as smooth as possible between adjacent pixels in the same frame image. The smaller the smoothness loss, the smoother the change of the optical flow field.

[0050] The gradient consistency loss L grad is defined as follows:

[0051]

[0052] where N is the total number of image pixels, is the value of the i-th pixel in the reference frame, v pred is the estimated optical flow field, is the optical flow vector of the i-th pixel in the estimated optical flow field, including the horizontal component and the vertical component, and represent calculating the gradient of the optical flow field in the x-axis and y-axis directions; the smaller the gradient consistency loss, the more consistent the change of the optical flow field with the image edge.

[0053] The final loss function is obtained by weighting the above losses:

[0054] L total =δL photo +βL smooth +γL grad

[0055] where δ, β, and γ are the hyperparameter weights of each loss term, controlling the contribution degree of each loss term in the total loss.

[0056] The present invention provides a cardiac motion feature extraction system based on optical flow estimation, including:

[0057] A preprocessing module for normalizing and cropping cardiac magnetic resonance data;

[0058] A feature extraction module, including a feature encoder and a context encoder, for extracting the features of the reference frame and the moving frame;

[0059] A correlation volume calculation module for constructing a correlation pyramid and extracting matching features;

[0060] A motion feature iterative enhancement module, including deformable convolution and a GMA (Global Motion Aggression) module, for refining the optical flow estimation;

[0061] A loss function optimization module for calculating the photometric consistency, smoothness, and gradient consistency losses and optimizing the model parameters.

[0062] Compared with the prior art, the advantages of the present invention are:

[0063] 1. The present invention adaptively captures local motion features at the edges of the cardiac structure through deformable convolution operations, effectively solving the problem of estimation errors of traditional optical flow methods in the edge region.

[0064] 2. The present invention introduces a Global Motion Aggression (GMA) module, which fuses global context features and local motion features through an attention mechanism to improve the optical flow estimation accuracy in large motion regions.

[0065] 3. The present invention designs a multi-loss function combining photometric consistency loss, smoothness loss, and gradient consistency loss to ensure the accuracy, smoothness, and edge consistency of the optical flow field, and improve the overall robustness of cardiac motion feature extraction.

[0066] 4. Based on an end-to-end framework of deep learning, the present invention does not require manual annotation of anatomical landmark points and can automatically extract high-precision motion features from cardiac magnetic resonance images, meeting the clinical requirements of high efficiency and accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 is the network structure proposed by the present invention.

[0068] Figure 2 is a schematic diagram of related volume calculation.

[0069] Figure 3 is the structure of the motion feature iterative enhancement module.

[0070] Figure 4 is the structure of the deformable convolution module.

[0071] Figure 5 is the structure of the Global Motion Aggression (GMA).

[0072] Figure 6 is the optical flow extraction result of the present invention. In Figure 6 , (a) is the reference frame image, (b) is the moving frame image, and (c) is the extracted optical flow arrow diagram. The direction of the arrow represents the cardiac motion direction, and the length represents the motion amplitude. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0073] The following embodiments will further illustrate the present invention with reference to the accompanying drawings. The embodiments of the present invention use the proposed optical flow estimation method to complete the extraction of cardiac motion features from cardiac magnetic resonance imaging data.

[0074] This embodiment adopts as Figure 1The network structure of the shown model divides two consecutive frames of images in time into a reference frame and a moving frame in sequence as input images and inputs them into the feature encoding module. The feature encoding module includes a feature encoder g θ (·) and a context encoder h λ (·), where both θ and λ are learnable network parameters; the extracted feature maps are input into the motion feature iterative enhancement module through the correlation volume calculation module, the attention module, and the initialized optical flow to optimize the optical flow estimation result. The specific structures of each module are described in detail in the specific implementation steps.

[0075] The embodiments of the present invention include the following steps:

[0076] 1) Preprocess the input four-dimensional spatio-temporal cardiac magnetic resonance imaging data, that is, scale the pixel values;

[0077] The four-dimensional spatio-temporal cardiac magnetic resonance data used is from the ACDC (Automated Cardiac Diagnosis Challenge) cardiac magnetic resonance dataset, with a total of 150 cases, 100 cases in the training set, and 50 cases in the test set; two consecutive frames of images in time are divided into a reference frame and a moving frame in sequence as input images; to reduce its impact on cardiac motion extraction and reduce the computational amount, it is necessary to preprocess the input four-dimensional cardiac magnetic resonance imaging data, that is, perform normalization processing and scale the pixel values to the range [0,1]; the process of normalization processing can be expressed as:

[0078]

[0079] where X is the input four-dimensional spatio-temporal cardiac magnetic resonance data, with a size of 256×216×1×10, A = 256 is the height of the image, B = 216 is the width of the image, S = 1 is the number of slices of the image, and fr = 10 is the number of time frames of the image; X norm is the data obtained after normalization processing, X min and X max are the minimum and maximum values of the pixels in the same image; traverse all time frames and the number of slices to crop the two-dimensional image with a size of A×B to obtain an image with a height and width of 128×128.

[0080] 2) Divide two consecutive frames of images in time into a reference frame and a moving frame, and input them into the feature encoding module for feature extraction. The feature encoding module includes a feature encoder and a context encoder; the feature encoder inputs the reference frame and the moving frame images, while the context encoder only inputs the reference frame image.

[0081] The feature encoder adopts a convolutional layer structure and encodes the input two frames of images into low-dimensional feature maps respectively; the encoder g θ(·) Encode the reference frame I1 and the moving frame I2 into features that are reduced by times in both length and width, namely g θ : where θ is a learnable network parameter. After the convolutional layer, there is a normalization layer and a rectified linear unit. The convolutional kernel size of the initial convolutional layer is 7×7; the feature encoder also includes 6 residual modules, each residual module contains two convolutional layers, and their convolutional kernel sizes are both 3×3. Each convolutional layer is followed by a normalization layer and a rectified linear unit; the feature encoding module also includes a context encoder h λ (·), λ is a learnable network parameter. The architecture of the context encoder is similar to that of the feature encoder. It extracts features from the reference frame I1, which is independent of the encoder g θ (·) respectively. The output feature extraction results are input into the correlation volume calculation module and the attention module respectively.

[0082] 3) Calculate the correlation between feature maps and construct a correlation pyramid through the correlation volume calculation schematic diagram shown as Figure 2 below, and extract feature maps from it to provide matching information for subsequent optical flow estimation; specifically:

[0083] Perform dot product on each corresponding pixel point of all the extracted feature pairs g θ (I1) and g θ (I2) respectively to obtain the correlation volume C. This process can be expressed as:

[0084]

[0085] where g θ (I1) and g θ (I2) are the extracted feature pairs, H = 128 and W = 128 are the height and width of the feature map respectively, and i, j, k, l, and h are the coefficients for traversing their respective dimensions.

[0086] The correlation volume calculation module consists of two parts: a correlation pyramid and a correlation lookup; in the correlation pyramid, use four pooling layers to pool the last two dimensions of the correlation volume C to construct a four-layer correlation pyramid {C 1 , C 2 , C 3 , C 4}; where the size of the Kth layer (K = 1, 2, 3, 4) of the correlation pyramid C K is In the correlation lookup operation L c L c extracts feature maps from the correlation pyramid; given the current optical flow estimate (f 1 , f 2), estimate the corresponding pixel U′=(u+f 1 (u),v+f 2 (v)); Define the search domain N(U′) of U′ r ;

[0087]

[0088] Where r is the search radius, dU is the offset, which belongs to the two-dimensional integer space, and ||·||1 is the L1 norm. For the relevant volume C k , in the search domain N(U′ / 2 K ) r Find the feature map (K=1,2,3,4) in the dataset, and finally concatenate the four feature maps into a correlation tensor.

[0089] 4) Use Figure 3 The motion feature iterative enhancement module shown in the figure continuously refines the optical flow estimation result through deformable convolution and global motion aggregation GMA (Global Motion Aggression) modules.

[0090] The motion feature iterative enhancement module inputs the initial estimation of optical flow and the context features processed by the attention module, and performs a deformable convolution operation on the two. The deformable convolution operation uses the following method: Figure 4 Operation shown; in a deformable convolution operation, a regular grid is sampled on the input feature map Increase the offset {Δp m |m=1,...,M}, where p m Represents traversal The position in the convolution kernel, M is the number of sampling points in the convolution kernel, that is, |·| is the number of elements in the set; for each position p0 on the output feature map Y, there is:

[0091]

[0092] in, For a regular network, Δp m is the offset, p m Represents traversal of a regular grid The position in the m ) is the convolution kernel at p m The value at , T is the input feature map.

[0093] The input feature map T is implemented by bilinear interpolation as follows:

[0094]

[0095] where p represents any position of the input feature map (p = p0 + p m + Δp m ), q traverses all integer spatial positions in the input feature map T, and G(·, ·) is a bilinear interpolation kernel.

[0096] The current optical flow information obtained after the deformable convolution operation forms a joint representation of the motion information and the spatial matching relationship with the correlation volume, and a global weighted modulation is imposed on this joint representation through the global motion aggregation GMA (Global Motion Aggression) as shown in Figure 5 In the global motion aggregation GMA (Global Motion Aggression) module, given the context feature and the motion feature where N = H × W and W are the height and width of the feature map respectively; represents the i-th context feature vector, represents the i-th motion feature vector, D c and D m are the numbers of the context feature vector and the motion feature vector respectively; the aggregated motion feature can be expressed as:

[0097]

[0098] where, is the updated i-th motion feature vector, represents the i-th motion feature vector, α is a learnable scalar parameter, φ and σ are the projection functions of the query vector, the key vector and the value vector respectively; f is a similarity attention function, which can be expressed as:

[0099]

[0100] where D is the number of channels of the feature map, φ and σ are the projection functions of the query vector, the key vector and the value vector respectively, and are expressed as:

[0101]

[0102] φ(x i ) = W key x i

[0103] σ(x i ) = W val x i

[0104] where is the weight matrix, Denote the i-th context feature vector.

[0105] The output of the Global Motion Aggression (GMA) module is the concatenation of the context feature, the motion feature, and the aggregated motion feature. The two-dimensional relative position embedding enables the attention map to depend not only on the self-similarity of the features but also on the relative position of the query points. After introducing the two-dimensional relative position embedding, the aggregated motion vector can be expressed as:

[0106]

[0107] Denote the i-th motion feature vector, α is a learnable scalar parameter, f is the similarity attention function, φ and σ are the projection functions of the query vector, the key vector, and the value vector respectively, and z j-i denotes the relative position embedding vector between positions j and i; the model learns the vertical and horizontal offset embedding vectors respectively and adds them to obtain z j-i .

[0108] Since the resolution of the predicted optical flow is that of the cropped image so upsample the optical flow to make its resolution the same as that of the cropped image.

[0109] 5) Optimize the model parameters based on the weighted sum of the photometric consistency loss, the smoothness loss, and the gradient consistency loss.

[0110] The loss function proposed by the present invention is as follows. For a pair of temporally consecutive image reference frame I1 and moving frame I2, the optical flow field v is used to inversely warp I1 to obtain the image I1'; the photometric consistency loss L photo is defined as follows:

[0111]

[0112] where, is the value of the i-th pixel in the estimated image obtained by inverse warping, is the value of the i-th pixel in the moving frame, N is the total number of image pixels, is the square of the L2 norm; the photometric consistency loss is used to measure the difference between the optical flow field warped image and the real image. The smaller the photometric consistency loss, the more accurate the optical flow estimation.

[0113] The smoothness loss L smooth is defined as follows:

[0114]

[0115] Where N is the total number of image pixels, v pred is the estimated optical flow field, is the optical flow vector of the i-th pixel of the estimated optical flow field, including the horizontal component and the vertical component; and represent taking the gradient of the optical flow field in the x-axis and y-axis directions; the smoothness loss is used to limit the discontinuity in the optical flow field, making the change of the optical flow as smooth as possible between adjacent pixels in the same frame image. The smaller the smoothness loss, the smoother the change of the optical flow field.

[0116] The gradient consistency loss L grad is defined as follows:

[0117]

[0118] Where N is the total number of image pixels, is the value of the i-th pixel of the reference frame, v pred is the estimated optical flow field, is the optical flow vector of the i-th pixel of the estimated optical flow field, including the horizontal component and the vertical component, and represent taking the gradient of the optical flow field in the x-axis and y-axis directions; the smaller the gradient consistency loss, the more consistent the change of the optical flow field with the edges of the image.

[0119] The final loss function is obtained by weighting the above losses:

[0120] L total = δL photo + βL smooth + γL grad

[0121] Where δ = 1, β = 0.05, and γ = 0.1 are the hyperparameter weights of each loss term, respectively, used to control the contribution degree of each loss term in the total loss.

[0122] The present invention is compared with several existing optical flow estimation methods (UnFlow, PWCNet, LiteFlowNet, and SpyNet) to evaluate the effectiveness of this method. The evaluation metrics include the mean squared error MSE (Mean Squared Error), peak signal-to-noise ratio PSNR (Peak Signal-to-Noise Ratio), and gradient error between the reference frame and the warped reference frame; the comparison of the performance metrics of the present invention and several optical flow estimation methods is shown in Table 1.

[0123] Table 1

[0124]

[0125] Note: Bold indicates the optimal method, and underlining indicates the sub-optimal method.

[0126] As can be seen from Table 1, the mean square error of the method of the present invention is reduced by 11.0% compared to the sub-optimal model SpyNet, indicating that the method of the present invention is more accurate in estimating the cardiac pixel displacement; the average peak signal-to-noise ratio of the method of the present invention reaches 29.72, which is 2.06% higher than that of the sub-optimal SpyNet and is the highest value among all the comparison models, indicating that the method of the present invention has a stronger signal recovery ability in the cardiac structure region, can more effectively suppress noise interference, and improve the image restoration quality; the average gradient error of the method of the present invention is 7.29, which is reduced by 9.57% compared to PWC-Net, indicating that the method of the present invention can more accurately predict the displacement field in the boundary regions such as the ventricular wall and valves.

[0127] Figure 6 is the optical flow extraction result of this method; in Figure 6 , (a) is the reference frame image, (b) is the moving frame image, and (c) is the extracted optical flow arrow diagram. The direction of the arrow represents the cardiac motion direction, and the length represents the motion amplitude; from Figure 6 it can be seen that the method proposed by the present invention can accurately capture the motion direction and amplitude, especially in the cardiac valve and myocardial regions, and its optical flow estimation result is more in line with the actual motion pattern.

[0128] The above embodiments are only preferred embodiments of the present invention and should not be considered as limiting the scope of implementation of the present invention. All equivalent changes and improvements made in accordance with the scope of the application of the present invention should still fall within the scope covered by the patent of the present invention.

Claims

1. A method for extracting cardiac motion features based on optical flow estimation, characterized in that Including the following steps: 1) Preprocess the input four-dimensional spatio-temporal cardiac magnetic resonance imaging data and scale the pixel values; 2) Divide two consecutive frames of images into a reference frame and a moving frame, and input them into the feature encoding module for feature extraction. The feature encoding module includes a feature encoder and a context encoder. The feature encoder encodes both frames of images, and the context encoder only encodes the reference frame; 3) Calculate the correlation between feature maps through the correlation volume calculation module, construct a correlation pyramid and extract feature maps from it to provide matching information for subsequent optical flow estimation; 4) The motion feature iterative enhancement module iteratively refines the optical flow estimation result through deformable convolution and the global motion aggregation GMA module; 5) Optimize the model parameters based on the weighted sum of the photometric consistency loss, smoothness loss, and gradient consistency loss.

2. The method for extracting cardiac motion features based on optical flow estimation according to claim 1, characterized in that In step 1), the preprocessing of the input four-dimensional spatio-temporal cardiac magnetic resonance imaging data, that is, normalization processing, scales the pixel values. The process of normalization processing is expressed as: Among them, X is the input four-dimensional spatio-temporal cardiac magnetic resonance data with a size of A×B×S×fr, where A is the height of the image, B is the width of the image, S is the number of image slices, and fr is the number of image time frames; X norm is the data obtained after normalization, X min and X max are the minimum and maximum values of the pixels in the same image; a two-dimensional image with a size of A×B is cropped to obtain an image with a height and width of H×W.

3. The method for extracting cardiac motion features based on optical flow estimation according to claim 1, wherein In step 2), the specific steps of the feature extraction are: The feature encoder adopts a convolutional layer structure to encode the two input frames of images into low-dimensional feature maps respectively; encoder g θ (·) encodes the reference frame I1 and the moving frame I2 into features that are respectively scaled down by a factor of where θ is a learnable network parameter, H and W are the height and width of the feature map respectively, ch is the number of channels, and the convolutional layer is followed by a normalization layer and a rectified linear unit; the feature encoder also includes J residual modules, each residual module contains two convolutional layers, and each convolutional layer is followed by a normalization layer and a rectified linear unit; the feature encoding module also includes a context encoder h λ (·), λ is a learnable network parameter, the architecture of the context encoder is similar to that of the feature encoder, extracts features from the reference frame I1, is independent of the encoder g θ (·) respectively, and the output feature extraction results are input into the correlation volume calculation module and the attention module.

4. The method for extracting cardiac motion features based on optical flow estimation according to claim 1, wherein In step 3), the specific steps of calculating the correlation between feature maps through the correlation volume calculation module, constructing a correlation pyramid and extracting feature maps from it to provide matching information for subsequent optical flow estimation are: For all pairs of the extracted features g θ (I1) and g θ (I2), the dot product is performed on each corresponding pixel point to obtain the correlation volume C, and this process is expressed as: Among them, g θ (I1) and g θ (I2) are the extracted feature pairs, H and W are the height and width of the feature map respectively, and i, j, k, l, and h are the coefficients for traversing their respective dimensions; The related volume calculation module consists of two parts: a correlation pyramid and a correlation lookup. In the correlation pyramid, the last two dimensions of the correlation volume C are pooled using four pooling layers to construct a four-layer correlation pyramid {C 1 , C 2 , C 3 , C 4}; where the size of the K-th layer of the correlation pyramid C K is In the correlation lookup operation L c , L c extracts feature maps from the correlation pyramid, where K = 1, 2, 3, 4. Given the current optical flow estimate (f 1 , f 2 ), the corresponding pixel point U' = (u + f 1 (u), v + f 2 (v)) of each pixel point U = (u, v) in the reference frame is estimated in the moving frame. Define the search domain N(U') of U' r where r is the search radius, dU is the offset, belonging to a two-dimensional integer space, ||·||1 is the L1 norm; for the correlation volume C k , Search for the feature map in the search domain N(U′ / 2 K ) r where K = 1, 2, 3, 4, and finally splice the four feature maps into a correlation tensor.

5. The method for extracting cardiac motion features based on optical flow estimation according to claim 1, wherein In step 4), the specific steps of the motion feature iterative enhancement module iteratively refines the optical flow estimation result through deformable convolution and the global motion aggregation GMA module are; The motion feature iterative enhancement module takes the initial optical flow estimate and the context features processed by the attention module as inputs, and performs deformable convolution operations on the two; in the deformable convolution operation, a regular grid sampled on the input feature map adds offsets {Δp m | m = 1,..., M}, where p m represents traversing the positions in, and M is the number of sampling points in the convolution kernel, that is |·| is the number of elements in the set; for each position p0 on the output feature map Y, there is: Among them, is a regular network, Δp m is the offset, p m represents traversing the positions in the regular grid in, w(p m ) is the value of the convolution kernel at p m , and T is the input feature map; The input feature map T is achieved by bilinear interpolation as follows: where \(p\) represents an arbitrary position in the input feature map, \(p = p_0 + p\) m +\Delta p m ; \(q\) traverses all integer spatial positions in the input feature map \(T\), and \(G(\cdot,\cdot)\) is a bilinear interpolation kernel; The current optical flow information obtained after the deformable convolution operation forms a joint representation of the motion information and the spatial matching relationship with the correlation volume, and a global weighted modulation is applied to this joint representation through the Global Motion Aggregation (GMA) module; in the GMA module, given the context features and the motion features where N = H × W and W are the height and width of the feature map respectively; represents the i-th context feature vector, represents the i-th motion feature vector, D c and D m are the numbers of the context feature vector and the motion feature vector respectively; the aggregated motion feature is expressed as: Among them, is the updated i-th motion feature vector, represents the i-th motion feature vector, and α is a learnable scalar parameter. φ and σ are projection functions for query vector, key vector, and value vector respectively; f is a similarity attention function, expressed as: where D is the number of channels of the feature map, φ and σ are the projection functions of the query vector, key vector, and value vector, respectively, which are expressed as: φ(x i ) = W key x i σ(x i ) = W val x i Among them, is the weight matrix, represents the i-th context feature vector; The output of the global motion aggregation (GMA) module is the concatenation of the context feature, the motion feature, and the aggregated motion feature After introducing the two-dimensional relative position embedding, the aggregated motion vector is expressed as: denotes the i-th motion feature vector, α is a learnable scalar parameter, and f is a similarity attention function, φ and σ are the projection functions for the query vector, key vector, and value vector respectively, and z j-i denotes the relative position embedding vector for positions j and i; the model learns the vertical and horizontal offset embedding vectors separately and adds them to obtain z j-i .

6. The method for extracting cardiac motion features based on optical flow estimation according to claim 1, characterized in that In step 5), the specific loss function for optimizing the model parameters based on the weighted sum of the photometric consistency loss, smoothness loss, and gradient consistency loss is as follows: For a pair of temporally consecutive image reference frames I1 and moving frame I2, the optical flow field v is used to perform inverse warping on I1 to obtain the image I1'; the photometric consistency loss L photo is defined as follows: wherein, is the value of the i-th pixel in the estimated image obtained by reverse deformation, is the value of the i-th pixel in the moving frame, N is the total number of image pixels, is the square of the L2 norm; the photometric consistency loss is used to measure the difference between the deformed image of the optical flow field and the real image. The smaller the photometric consistency loss, the more accurate the optical flow estimation; Smoothness loss L smooth is defined as follows: where N is the total number of image pixels, and v pred is the estimated optical flow field, and is the optical flow vector of the i-th pixel of the estimated optical flow field, including the horizontal and vertical components; and represent the gradient calculation of the optical flow field in the x-axis and y-axis directions; the smoothness loss is used to limit the discontinuity in the optical flow field, making the change of the optical flow between adjacent pixels in the same frame of image as smooth as possible. The smaller the smoothness loss, the smoother the change of the optical flow field; Gradient consistency loss L grad is defined as follows: where N is the total number of image pixels, is the value of the i-th pixel in the reference frame, v pred is the estimated optical flow field, is the optical flow vector of the i-th pixel in the estimated optical flow field, including the horizontal component and the vertical component, and represent taking the gradient of the optical flow field in the x-axis and y-axis directions; the smaller the gradient consistency loss, the more consistent the change of the optical flow field is with the edges of the image; The final loss function is obtained by weighting the above losses: L total = δL photo + βL smooth + γL grad where δ, β, and γ are the hyperparameter weights of each loss term, respectively, controlling the contribution degree of each loss term in the total loss.

7. A cardiac motion feature extraction system based on optical flow estimation, characterized in that Including: A preprocessing module for normalizing and cropping cardiac magnetic resonance data; A feature extraction module, including a feature encoder and a context encoder, for extracting features of the reference frame and the moving frame; A correlation volume calculation module for constructing a correlation pyramid and extracting matching features; A motion feature iterative enhancement module, including a deformable convolution and a GMA module, for refining optical flow estimation; A loss function optimization module for calculating the photometric consistency, smoothness, and gradient consistency losses and optimizing the model parameters.

Images