Motion vector estimation method and system for ultrasonic imaging, medium and equipment

Abstract

The invention discloses a motion vector estimation method and system for ultrasonic imaging, a medium and equipment. The method comprises the following steps: sequentially carrying out time high-pass filtering and spatial difference Gaussian band-pass filtering on an initial signal; constructing Reichardt correlators in a plurality of directions in the IQ complex domain, and obtaining a plurality of correlation energies in the plurality of directions; a multi-scale parameter set is set, the optimal speed of each scale is obtained based on the pixel physical spacing and the frame rate, the displacement under each scale is mapped into the optimal speed, fusion is carried out with the related energy as the weight, and omnidirectional speed estimation is output and comprises the speed amplitude and the direction unit vector; and obtaining an estimated axial speed of the initial signal based on a Kasai autocorrelation method, carrying out adaptive weighted fusion on the omnidirectional speed estimation and the axial speed, and outputting a motion vector and confidence. According to the invention, the estimation performance of the motion vector is improved, and the problem of insensitivity to low-speed weak current signals in the prior art is solved.

Description

Technical Field

[0001] This invention relates to the field of ultrasound imaging and signal processing technology, and specifically to a method, system, medium, and device for estimating motion vectors in ultrasound imaging. Background Technology

[0002] In the field of medical ultrasound imaging, motion vector estimation is crucial for the dynamic analysis of tissues and blood flow. Currently, the mainstream clinical velocity estimation techniques mainly include the following categories: First, Doppler techniques, such as color Doppler and pulsed Doppler, which estimate the velocity component along the ultrasound beam axis by measuring the phase change of the echo signal; Second, video processing methods based on B-mode images, such as block matching or optical flow methods, which infer motion by analyzing the grayscale changes of pixels or regions between consecutive frames.

[0003] In existing technologies, Doppler technology, limited by its physical principles, is highly sensitive to the angle between the motion direction and the beam direction, leading to a severe underestimation of non-axial motion components. Furthermore, its pulse repetition frequency (PRF) and Nyquist limit restrict its velocity upper limit to PRF / Nyquist. Insensitivity to low-speed, weak flow is primarily due to wall filtering and low SNR. On the other hand, B-mode video-based methods are susceptible to pseudo-scintillation interference from the inherent speckle noise of ultrasound, easily generating spurious motion vectors (pseudo-flows) under low signal-to-noise ratio conditions. While multi-angle emission or vector flow imaging (VFI) techniques can mitigate the angle-dependent problem to some extent, these approaches typically require significant modifications to ultrasound equipment and beam sequences, placing high demands on system computing resources and making it difficult to achieve low-cost, high-frame-rate real-time imaging in clinical settings.

[0004] Therefore, a new technological solution is needed. Summary of the Invention

[0005] In view of this, embodiments of the present invention provide a method, system, medium, and device for estimating motion vectors in ultrasound imaging, to at least solve the problem of insensitivity to low-speed weak flow signals in the prior art.

[0006] The embodiments of the present invention provide the following technical solutions: This invention provides a motion vector estimation method for ultrasound imaging, comprising: Acquire ultrasonic echo complex signals or their envelope sequences arranged in a time frame sequence as the initial signal; The initial signal is subjected to time high-pass filtering and spatial difference Gaussian band-pass filtering in sequence to suppress static background and enhance the spatial contrast of moving targets; Multiple Reichardt correlators in multiple directions are constructed in the IQ complex domain. Each Reichardt correlator is subjected to signal delay, conjugate multiplication, and antisymmetric differential operation to obtain multiple correlation energies in the multiple directions. A multi-scale parameter group is set, and the optimal speed for each scale is obtained based on the pixel physical spacing and frame rate. The displacement under each scale is mapped to the optimal speed, and the corresponding optimal speed amplitude is fused with the relevant energy as weight to output an omnidirectional speed estimate. The omnidirectional speed estimate includes speed amplitude and direction unit vector. The estimated axial velocity of the initial signal is obtained based on the Kasai autocorrelation method. The omnidirectional velocity estimate is adaptively weighted and fused with the axial velocity, and the motion vector and confidence level are output. The fusion weight is determined based on the Kasai autocorrelation amplitude and the correlation energy.

[0007] Furthermore, after performing the time high-pass filtering on the initial signal, an adaptive threshold is determined based on the interquartile range of the local intensity change and the median of the local intensity of the filtered initial signal. The inter-frame intensity difference of the initial signal is decomposed into a signal enhancement channel and a signal attenuation channel, and a dual-channel gating weight is constructed to suppress pseudo-intensity changes caused by coherent speckle. The adaptive threshold τ = κ·IQR_local + β·median, where κ and β are adjustment coefficients, and their values ​​range from κ∈[0.02, 0.15] and β∈[0, 0.2], respectively.

[0008] Furthermore, the time-pass filtering is implemented using the formula Y(t)=X(t)-αX(t-1); Where X(t) is the current frame signal, X(t-1) is the previous frame signal, and α is the filter coefficient, which takes values ​​in the range of [0.85, 0.98].

[0009] Furthermore, the spatial difference Gaussian bandpass filtering is achieved by the difference between two Gaussian kernels with standard deviations of σ1 and σ2, where σ1 ∈ [0.4, 0.8] pixels, σ2 ∈ [1.0, 1.8] pixels, and σ1 < σ2.

[0010] Furthermore, the spatial difference Gaussian bandpass filtering is implemented using a split convolution method, which decomposes the two-dimensional Gaussian convolution into sequentially executed horizontal and vertical one-dimensional convolutions.

[0011] Furthermore, the multiple directions include four directions: left, right, up, and down; wherein, the rightward related energy is calculated using the following formula: R→=Re{X(x,y,t)·X*(x+Δx,y,t+Δt)-X(x+Δx,y,t)·X*(x,y,t+Δt)}, Where X() is the signal, X*() is its conjugate, and Re represents taking the real part; the correlation energy in the left, up, and down directions is calculated by symmetrically adjusting the displacement parameters (Δx, Δy).

[0012] Furthermore, the multi-scale parameter set includes the number of horizontal displacement pixels Δx, the number of vertical displacement pixels Δy, and the number of time delay frames Δt; the preferred velocity v_pref is calculated using the following formula: ; Where Δ_x and Δ_y are the physical pixel spacing in the horizontal and vertical directions, respectively, and fps is the imaging frame rate.

[0013] Furthermore, before adaptively weighting and fusing the omnidirectional velocity estimate and the axial velocity, the two are first aligned in terms of time length. The alignment time window length is T_align=min(T_corr, T_dop), where T_corr is the time length used to calculate the omnidirectional velocity estimate and T_dop is the time length used to calculate the axial velocity.

[0014] Furthermore, in the adaptive weighted fusion, the fusion weight w_d assigned to the axial velocity is determined by the following function: w_d=γ·C_d / (γ·C_d+(1-γ)·C_c); Where C_d is the normalized Kasai autocorrelation magnitude, C_c is the normalized correlation energy, and γ is an adjustable parameter with a value range of [0.3, 0.8].

[0015] This invention provides a motion vector estimation system for ultrasound imaging, comprising: The time-domain gating module is used to acquire ultrasonic echo complex signals or their envelope sequences arranged in a time frame sequence as initial signals; and to sequentially perform time high-pass filtering on the initial signals to suppress static background. The spatial bandpass module is used to perform spatial difference Gaussian bandpass filtering on the time-domain processed signal received from the time-domain gating module to enhance the spatial contrast of the moving target. An IQ domain directional correlator array is used to construct multiple directional Reichardt correlators in the IQ complex domain. Each Reichardt correlator is processed by signal delay, conjugate multiplication, and antisymmetric differential operation to obtain multiple correlation energies in the multiple directions. The multi-scale velocity calibration and fusion module is used to receive the correlation energy output by the IQ domain direction correlator array, set a multi-scale parameter group, and obtain the preferred velocity for each scale based on the pixel physical spacing and frame rate. The displacement under each scale is mapped to the preferred velocity, and the corresponding preferred velocity amplitude is fused with the correlation energy as weight to output an omnidirectional velocity estimate. The omnidirectional velocity estimate includes velocity amplitude and direction unit vector. The Kasai Doppler module is used to receive the initial signal, calculate and output the estimated axial velocity using the Kasai autocorrelation method; The adaptive fusion and confidence module is used to receive the estimated axial velocity and perform adaptive weighted fusion of the omnidirectional velocity estimate and the axial velocity, wherein the fusion weight is determined based on the Kasai autocorrelation amplitude and the correlation energy, and finally outputs the motion vector and confidence.

[0016] Furthermore, the spatial bandpass module employs a split convolution method to achieve spatial difference Gaussian bandpass filtering; or The IQ domain directional correlator array consists of programmable delay line units and complex multiplier / accumulator units. The programmable delay line units perform signal delay operations, and the complex multiplier / accumulator units perform numerical calculations of conjugate multiplication and antisymmetric difference.

[0017] The present invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, performs the steps of any of the motion vector estimation methods described above.

[0018] The present invention also provides an electronic device, comprising: a processor, a memory, and a bus, wherein the memory stores machine-readable instructions executable by the processor, and when the electronic device is running, the processor communicates with the memory via the bus, and the machine-readable instructions are executed by the processor to perform the steps of any of the motion vector estimation methods described above.

[0019] Compared with the prior art, the beneficial effects that the at least one technical solution adopted in the embodiments of the present invention can achieve include at least: This invention improves the estimation performance of motion vectors through a collaborative technology chain of spatiotemporal preprocessing, omnidirectional velocity extraction, and axial velocity fusion; enhances the sensitivity to weak and slow currents through temporal high-pass and spatial DoG band-pass filtering; furthermore, it improves non-axial detection capability through an IQ complex domain four-directional Reichardt correlator, and fully captures lateral motion components using complex signal phase information. Attached Figure Description

[0020] To more clearly illustrate the technical solutions of the embodiments of this application, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1-2 This is a flowchart of a motion vector estimation method for ultrasound imaging according to an embodiment of the present invention; Figure 3 This is a schematic diagram of the timing of time-pass filtering and IQR-based ON / OFF gating in an embodiment of the present invention; Figure 4 This is a schematic diagram of the IQ domain four-directional Reichardt correlator structure and signal path according to an embodiment of the present invention; Figure 5 This is a schematic diagram of the adaptive fusion weighting function of Kasai autocorrelation and correlator energy in an embodiment of the present invention; Figure 6 This is a flowchart illustrating the ON / OFF gate control process according to an embodiment of the present invention. Figure 7 This is a speed / direction overlay diagram and confidence diagram of a software embodiment of the present invention; Figure 8 The ROC / AUC and minimum detectable speed (MDS) curves are shown in the comparison and ablation experiments of this invention. Detailed Implementation

[0022] The embodiments of this application will now be described in detail with reference to the accompanying drawings.

[0023] The following specific examples illustrate the implementation of this application. Those skilled in the art can easily understand other advantages and effects of this application from the content disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. This application can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of this application. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0024] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this application, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number and aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0025] It should also be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of this application. The drawings only show the components related to this application and are not drawn according to the actual number, shape and size of the components in the actual implementation. In the actual implementation, the form, quantity and proportion of each component can be arbitrarily changed, and the layout of the components may also be more complex.

[0026] Additionally, specific details are provided in the following description to facilitate a thorough understanding of the examples. However, those skilled in the art will understand that practice can be carried out without these specific details.

[0027] Color / pulse Doppler is commonly used for velocity estimation, but it is limited by axial angle and PRF / Nyquist, and is insensitive to non-axial components and weak flows at low speeds (e.g., <5 cm / s). B-mode video methods based on block matching or optical flow are susceptible to speckle "flickering" and low signal-to-noise ratio, resulting in spurious flows. While multi-angle / vector flow imaging (VFI) can improve the angle problem, it requires significant resources, equipment modifications, and sophisticated beamforming schemes, leading to high costs for real-time clinical implementation. Therefore, there is an urgent need for a resource-efficient solution that can robustly output motion direction and velocity vectors under single-angle or limited-angle, low SNR conditions. Based on this, the embodiments of this specification propose a processing solution: such as Figure 1As shown, this invention first performs a time high-pass on the input ultrasonic frame sequence and constructs an ON / OFF dual-channel gating based on the local interquartile range (IQR) to suppress speckle. Then, a spatial bandpass is achieved using the difference-of-Gaussian (DoG) method. Subsequently, a four-directional Reichardt correlator is constructed in the IQ complex domain and multi-scale processing is performed to map (Δx, Δy, Δt) to a preferred velocity and fuse it using energy weighting to obtain the velocity and direction. Simultaneously, the Kasai autocorrelation method is used to estimate the axial velocity, and adaptive weighted fusion is performed based on a function of the autocorrelation amplitude and correlation energy to output the motion vector and confidence level. This scheme exhibits higher sensitivity and robustness to slow flow, weak flow, and non-axial motion under single-angle or finite-angle, low signal-to-noise ratio conditions. Furthermore, it employs a hardware-friendly implementation using separable convolution, delay lines, and complex multiplication / accumulation, making it suitable for real-time deployment on FPGA / SoC / GPU.

[0028] Specifically, the motion vector estimation method of the present invention can be summarized as a cat's-eye spatiotemporal signal chain: Time-based high-pass + ON / OFF gating: Perform a first-order differential high-pass Y(t)=X(t)-αX(t-1) on the frame sequence, and decompose the intensity change into ON (enhancement) and OFF (weakening) channels based on the adaptive threshold of local IQR to form gating weights to suppress speckle flicker; Spatial bandpass (DoG): Enhances the spatial contrast between moving edges and speckle clusters; IQ domain four-directional Reichardt correlator (multi-scale): Construct four-directional correlation units (→, ←, ↑, ↓) in the complex IQ domain using "delay-conjugate multiplication-antisymmetric difference" to perform directional energy response on the (Δx, Δy, Δt) multi-scale group; Velocity physics calibration and fusion: The displacement / delay at each scale is mapped to the preferred velocity v_pref, and the velocity and direction are fused with relevant energy weights; Kasai Doppler (axial) × correlator (transverse / full field) adaptive fusion: Confidence is constructed using the first-order autocorrelation amplitude |R1| and correlation energy, and then adaptively fused according to the weight function; Time length alignment and post-processing: T_align=min(T_corr,T_dop), outputs velocity vector, direction angle and confidence level.

[0029] The technical solutions provided by the various embodiments of this application are described below with reference to the accompanying drawings.

[0030] like Figure 1-2 As shown, the present invention provides a motion vector estimation method for ultrasound imaging, comprising: Step S102: Acquire the ultrasonic echo complex signal or its envelope sequence arranged in time frame sequence as the initial signal.

[0031] The ultrasound echo complex signal (IQ signal) is a complex signal consisting of a real part (I) and an imaginary part (Q). It can be represented as A·e^(jφ), where A is the amplitude and φ is the phase. In ultrasound systems, it is obtained by quadrature demodulation of the radio frequency (RF) signal and is standard data for advanced hemodynamic analysis and vector velocity estimation. The envelope sequence is the result of converting the complex form of the IQ signal into a real form, which is the B-mode (luminance mode) grayscale image we usually see.

[0032] Among them, an ultrasound probe can be used to emit ultrasound waves into the human body and receive the echo signals scattered or reflected back from tissues and blood flow.

[0033] The initial signals are arranged frame by frame in chronological order to form a three-dimensional dataset.

[0034] Step S102 is used to provide the raw inputs required by all subsequent processing modules (time-domain filtering, spatial filtering, correlator, Doppler estimation).

[0035] Step S104: Perform time high-pass filtering and spatial difference Gaussian band-pass filtering on the initial signal in sequence to suppress static background and enhance the spatial contrast of moving targets.

[0036] Combination Figure 3 As shown, the initial signal needs to be subjected to time high-pass filtering first, and then the pseudo signal needs to be processed before spatial difference Gaussian band-pass filtering is performed.

[0037] The time-pass filtering is achieved by the formula Y(t)=X(t)-αX(t-1); where X(t) is the current frame signal, X(t-1) is the previous frame signal, and α is the filtering coefficient, which takes values ​​in the range of [0.85, 0.98].

[0038] Specifically, because the background organization is static or slowly changing, its signals in adjacent frames are very similar, and the result after subtraction is close to zero, thus filtering out frames with static backgrounds using the above formula; because the moving target changes its position or shape between frames, the subtraction will produce a significant differential signal, thus highlighting it.

[0039] Therefore, time-pass filtering is used to filter out static background signals that do not change over time, providing a "clean" time-series signal containing only motion information for subsequent motion estimation steps.

[0040] After performing time-pass filtering on the initial signal, the following steps are also included: Step S104a: Based on the interquartile range and median local intensity variation of the filtered initial signal, an adaptive threshold is determined, and the inter-frame intensity difference of the initial signal is decomposed into a signal enhancement channel and a signal attenuation channel. A dual-channel gating weight is constructed to suppress pseudo-intensity variations caused by coherent speckles. The adaptive threshold is as follows: τ=κ·IQR_local+β·median, κ∈[0.02, 0.15], β∈[0, 0.2]; Wherein, κ and β are adjustment coefficients, with values ​​ranging from κ ∈ [0.02, 0.15] to β ∈ [0, 0.2].

[0041] Among them, the interquartile range of local intensity variation is used to measure the dispersion of signal changes in a local area; the median local intensity represents the baseline level of signal intensity in that area; and the adaptive threshold is used to provide a dynamic judgment criterion related to the statistical characteristics of the local image, avoiding the problem of using a fixed threshold that is too sensitive in some areas (introducing noise) or too insensitive in other areas (losing the true signal).

[0042] The signal enhancement channel retains only the portion of the inter-frame intensity difference that is greater than zero, while the signal attenuation channel retains only the portion of the inter-frame intensity difference that is less than zero. The comparison result between the inter-frame intensity difference and the adaptive threshold generates a gating weight for each channel. Only when the absolute value of the inter-frame intensity difference is greater than the adaptive threshold will it be assigned a higher weight.

[0043] The gating weight is used to regulate the "trust level" of subsequent processing steps towards the signal of the current pixel. A weight of 1 indicates that the change at this point is completely trusted as real motion; a weight of 0 indicates that the change is completely ignored and considered noise; and a weight in between indicates partial trust. This weight is directly multiplied by the signal at that point, thereby suppressing noise at the source.

[0044] Step S104a is located after time high-pass filtering and is used to distinguish between real motion signals and pseudo intensity changes caused by coherent speckle noise. It processes positive and negative signal changes separately through dual gates, enabling the system to process complex motion patterns more precisely.

[0045] Spatial difference Gaussian bandpass filtering is achieved by the difference between two Gaussian kernels with standard deviations σ1 and σ2, where σ1 ∈ [0.4, 0.8] pixels, σ2 ∈ [1.0, 1.8] pixels, and σ1 < σ2.

[0046] Among them, the spatial difference Gaussian bandpass filter can form a DoG filter by calculating the difference between two Gaussian kernels with different standard deviations (σ1 and σ2, and σ1<σ2).

[0047] Specifically, a small standard deviation Gaussian kernel (σ1) is used to preserve fine edges and details (which may contain noise), while a large standard deviation Gaussian kernel (σ2) is used to capture large-scale background structure. Subtracting the two cancels out the large-scale uniform background while preserving and enhancing features between the two scales.

[0048] Spatial difference Gaussian bandpass filtering weakens uniform regions in the spatial dimension and sharpens the boundaries and internal texture structure of moving targets, thereby improving the spatial recognizability of moving targets.

[0049] Furthermore, the spatial difference Gaussian bandpass filter is implemented using a split convolution method, which decomposes the two-dimensional Gaussian convolution into sequentially executed horizontal and vertical one-dimensional convolutions.

[0050] Step S104 achieves coordinated preprocessing of "time domain first, then spatial domain" by sequentially performing time high-pass filtering and spatial difference Gaussian band-pass filtering. The time filtering first "liberates" the signal from the static background, and then the spatial filtering further enhances the spatial characteristics of these moving targets, thereby significantly improving the contrast between the moving targets and the background in the output signal. This provides a high-quality, high signal-to-noise ratio input for the subsequent IQ domain directional correlator, directly improving the system's ability to detect and estimate weak, non-axial motion.

[0051] Step S106: Construct multiple Reichardt correlators in the IQ complex domain. Each Reichardt correlator is processed by signal delay, conjugate multiplication, and antisymmetric differential operation to obtain multiple correlation energies in multiple directions.

[0052] Combination Figure 4 As shown, the Reichardt correlator is a motion detection model that primarily senses motion in a specific direction by comparing the signals from two adjacent spatial points at different times. Signal delay involves applying a small time delay Δt to the signal from one spatial point, simulating the time required for a motion signal to propagate from one point to another. Conjugate multiplication involves performing a complex multiplication of the undelayed signal with the conjugate of the delayed signal. In signal processing, this is equivalent to calculating the cross-correlation of the two signals, the result of which is a complex number whose phase contains information about the phase difference between the two signals, and the phase difference is extremely sensitive to motion. Antisymmetric difference calculates the difference between the correlation results in two opposite directions (e.g., right-hand correlation minus left-hand correlation). This is to impart directional selectivity to the response.

[0053] The multiple directions include four directions: left, right, up, and down; the rightward energy is calculated using the following formula: R→=Re{X(x,y,t)·X*(x+Δx,y,t+Δt)-X(x+Δx,y,t)·X*(x,y,t+Δt)}, Where X() is the signal, X*() is its conjugate, and Re represents taking the real part; the correlation energy in the left, up, and down directions is calculated by symmetrically adjusting the displacement parameters (Δx, Δy).

[0054] Step S106 is used to detect and quantify the motion intensity (correlation energy) in the four basic directions of up, down, left, and right. This allows the Reichardt correlators in the four directions to sense and measure the lateral motion component perpendicular to the ultrasonic beam direction, overcoming the fundamental limitation of traditional Doppler technology, which can only measure axial velocity. Furthermore, step S106 outputs a clear "correlation energy" value for each direction. The magnitude of this value represents the intensity of the motion signal in that direction, and its sign indicates the direction of motion. This lays the foundation for subsequent fusion to generate accurate vector directions.

[0055] Step S106, by operating in the IQ complex domain, can utilize phase information, making it highly sensitive to extremely slow motion (slow flow) and real motion under low signal-to-noise ratio conditions, with performance far superior to traditional optical flow methods operating on B-mode images.

[0056] Step S108: Set a multi-scale parameter group, and obtain the preferred speed for each scale based on the pixel physical spacing and frame rate. Map the displacement under each scale to the preferred speed, and fuse the corresponding preferred speed amplitude with the relevant energy as weight to output the omnidirectional speed estimate. The omnidirectional speed estimate includes speed amplitude and direction unit vector.

[0057] The multi-scale parameter set includes the number of horizontal displacement pixels Δx, the number of vertical displacement pixels Δy, and the number of time delay frames Δt; the preferred velocity v_pref is calculated using the following formula: ; Where Δ_x and Δ_y are the physical pixel spacing in the horizontal and vertical directions, respectively, and fps is the imaging frame rate.

[0058] The multi-scale parameter set (Δx, Δy, Δt) is a series of pre-defined combinations covering different displacement ranges and time delays. For example, one set of parameters may be for fast flow with large displacements, while another set may be for slow flow with small displacements.

[0059] Specifically, a multi-scale parameter set is set (including the number of horizontal displacement pixels Δx, the number of vertical displacement pixels Δy, and the number of time delay frames Δt). Based on the pixel physical spacing (horizontal direction Δ_x, vertical direction Δ_y) and the frame rate (fps), the displacement at each scale is mapped to the optimal velocity for the corresponding scale. The optimal velocity amplitude for each scale is determined by the formula... The calculated direction vector is determined by normalizing the displacement direction (Δx, Δy) at that scale. Then, using the relevant energy corresponding to each scale as the weight, the preferred velocity amplitudes of all scales are weighted and summed (i.e., final velocity amplitude = Σ(relevant energy_i × v_pref_i) / Σrelevant energy_i), and the direction unit vectors of all scales are weighted and summed again and normalized (i.e., final direction unit vector = [Σ(relevant energy_i × direction vector_i)] / Σrelevant energy_i). The omnidirectional velocity estimate is output through the above fusion, and the omnidirectional velocity estimate includes the velocity amplitude and the direction unit vector.

[0060] The omnidirectional velocity estimation is the final output of this step. It integrates information from all directions and all scales, representing the motion state of each pixel in a single vector form. This vector consists of two parts: a velocity magnitude (a scalar representing speed) and a direction unit vector (representing only direction, with a length of 1), which fully describes the velocity and direction of the motion.

[0061] Step S106 uses multi-scale detection to ensure that the system can simultaneously and sensitively detect both fast and slow flows. Furthermore, through energy-weighted fusion, it can automatically filter and fuse the most reliable and consistent final estimate from multiple trial velocities, thereby improving the accuracy and robustness of the estimate.

[0062] Step S110: Obtain the estimated axial velocity of the initial signal based on the Kasai autocorrelation method, perform adaptive weighted fusion of the omnidirectional velocity estimate and the axial velocity, and output the motion vector and confidence level. The fusion weight is determined based on the Kasai autocorrelation amplitude and correlation energy.

[0063] Among them, the Kasai autocorrelation method is a Doppler velocity estimation algorithm. It estimates the phase change rate by calculating the autocorrelation function between adjacent ultrasonic echo signals, thereby calculating the velocity component along the ultrasonic beam direction. Its advantages are high computational efficiency and accurate axial velocity measurement, making it a standard algorithm for color Doppler imaging.

[0064] Before adaptively weighting and fusing the omnidirectional velocity estimate and the axial velocity, the two are first aligned in terms of time length. The alignment time window length is T_align=min(T_corr, T_dop), where T_corr is the time length used to calculate the omnidirectional velocity estimate and T_dop is the time length used to calculate the axial velocity.

[0065] The time window length formula determines a common effective time window and achieves synchronization by truncating longer data sequences, preventing errors caused by time misalignment.

[0066] Furthermore, in adaptive weighted fusion, the fusion weight w_d assigned to the axial velocity is determined by the following function: w_d=γ·C_d / (γ·C_d+(1-γ)·C_c); Where C_d is the normalized Kasai autocorrelation magnitude, which is the confidence level of the axial velocity path; C_c is the normalized correlation energy, which is the confidence level of the omnidirectional velocity path; and γ is an adjustable parameter with a value range of [0.3, 0.8], used to identify which path the overall system is biased towards.

[0067] Combination Figure 5-8 As shown, step S110 ensures temporal consistency of the fusion through time alignment, which is used to guarantee the smoothness and accuracy of the final motion vector in time. The adaptive weight is used to dynamically determine which path the final result depends on more, so as to obtain the optimal estimate in different regions such as the center of the blood vessel (axial flow dominates) and the edge / bend (lateral flow is significant), thereby enabling the output of the optimal motion vector.

[0068] The motion vector estimation method of the present invention has the following advantages: it is sensitive to weak / slow flow, has strong non-axial capability, high robustness (IQ domain correlation / conjugate multiplication), resource-friendly (DoG-separated convolution, delay line + complex multiplication MAC pipeline), and has strong system-level availability (it can be deployed on existing systems via software / firmware upgrades).

[0069] This invention discloses a motion vector estimation system for ultrasound imaging that combines the axial accuracy of the Kasai method with the lateral sensitivity of the Reichardt correlator, and simultaneously outputs a confidence map that quantifies the reliability of the velocity estimate for each pixel. Doctors can use this system to determine which blood flow signals are robust and reliable, and which may have uncertainties, greatly improving the reliability of clinical diagnosis.

[0070] The present invention also provides a motion vector estimation system for ultrasound imaging, comprising a time-domain gating module, a spatial bandpass module, an IQ-domain direction correlator array, a multi-scale velocity calibration and fusion module connected in sequence, and a Kasai Doppler module connected in sequence to the multi-scale velocity calibration and fusion module and the adaptive fusion and confidence module, respectively, wherein the IQ-domain direction correlator array is also connected in sequence to the adaptive fusion and confidence module.

[0071] The time-domain gating module is used to acquire ultrasonic echo complex signals or their envelope sequences arranged in a time frame sequence as the initial signal; the initial signal is sequentially subjected to time high-pass filtering to suppress static background; and an adaptive threshold is determined based on the interquartile range of local intensity changes and the median of local intensity of the filtered signal. The inter-frame intensity difference is decomposed into a signal enhancement channel and a signal attenuation channel, and a dual-channel gating weight is constructed to suppress pseudo-intensity changes caused by coherent speckle. The time-domain processed signal is then output to the spatial bandpass module.

[0072] The spatial bandpass module is used to perform spatial difference Gaussian bandpass filtering on the time-domain processed signal output by the receiving time-domain gating module to enhance the spatial contrast of moving targets.

[0073] Specifically, the spatial bandpass module performs spatial difference Gaussian bandpass filtering by using the difference between two Gaussian kernels with standard deviations of σ1 (σ1∈[0.4, 0.8] pixels) and σ2 (σ2∈[1.0, 1.8] pixels, and σ1<σ2) to enhance the spatial contrast of the moving target and outputs the spatially processed signal to the IQ domain direction correlator array.

[0074] The IQ domain directional correlator array is used to construct multiple directional Reichardt correlators in the IQ complex domain. Each Reichardt correlator is obtained by signal delay, conjugate multiplication, and antisymmetric differential operation, resulting in multiple correlation energies in multiple directions.

[0075] Specifically, the IQ domain directional correlator array receives the spatially processed signal output by the spatial bandpass module, constructs Reichardt correlators in the left, right, up, and down directions in the IQ complex domain, and each Reichardt correlator sequentially performs signal delay, conjugate multiplication, and antisymmetric difference operations to calculate and output the correlation energy in the four directions, which is then sent to the multi-scale velocity calibration and fusion module and the adaptive fusion and confidence module, respectively.

[0076] The multi-scale velocity calibration and fusion module is used to receive the correlation energy output by the IQ domain direction correlator array, set the multi-scale parameter group, and obtain the optimal velocity for each scale based on the pixel physical spacing and frame rate. The displacement under each scale is mapped to the optimal velocity, and the correlation energy is used as the weight to perform fusion and output the omnidirectional velocity estimate. The omnidirectional velocity estimate includes the velocity amplitude and the direction unit vector.

[0077] Specifically, the multi-scale velocity calibration and fusion module receives the correlation energy output by the IQ domain direction correlator array, sets a multi-scale parameter set including the number of horizontal displacement pixels Δx, the number of vertical displacement pixels Δy, and the number of time delay frames Δt, calculates the preferred velocity for each scale according to the preferred velocity formula based on the physical pixel spacing (Δ_x, Δ_y) in the horizontal and vertical directions and the imaging frame rate (fps), fuses the preferred velocities of each scale with the correlation energy as weight, and outputs an omnidirectional velocity estimate containing velocity amplitude and direction unit vector to the adaptive fusion and confidence module. The Kasai Doppler module is used to receive the initial signal, calculate and output the estimated axial velocity using the Kasai autocorrelation method.

[0078] The adaptive fusion and confidence module is used to receive the estimated axial velocity and perform adaptive weighted fusion of the omnidirectional velocity estimate and the axial velocity. The fusion weight is determined based on the Kasai autocorrelation amplitude and correlation energy, and finally outputs the motion vector and confidence.

[0079] Specifically, the adaptive fusion and confidence module receives the omnidirectional velocity estimate and omnidirectional velocity calculation time T_corr output by the multi-scale velocity calibration and fusion module, the estimated axial velocity and axial velocity calculation time T_dop output by the Kasai Doppler module, and the correlation energy output by the IQ domain direction correlator array. It first aligns the omnidirectional velocity estimate and the estimated axial velocity by time length according to T_align=min(T_corr, T_dop), and then, based on the Kasai autocorrelation amplitude (normalized to C_d) and the correlation energy (normalized to C_c), it calculates the axial velocity estimate according to the formula... Determine the fusion weights, perform adaptive weighted fusion on the aligned omnidirectional velocity estimate and the estimated axial velocity, and finally output the motion vector and confidence level.

[0080] Furthermore, the spatial bandpass module employs a split convolution method to implement spatial difference Gaussian bandpass filtering, decomposing the two-dimensional Gaussian convolution into sequentially executed horizontal one-dimensional convolution and vertical one-dimensional convolution, reducing computational complexity to adapt to real-time processing.

[0081] Furthermore, the IQ domain directional correlator array consists of programmable delay line units and complex multiplier / accumulator units. The programmable delay line units perform signal delay operations, while the complex multiplier / accumulator units perform numerical calculations of conjugate multiplication and antisymmetric difference, adapting to the pipelined data processing of FPGA / SoC.

[0082] This application also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, can perform the above-described actions. Figure 1The steps of the motion vector estimation method in the illustrated embodiment can be found in the method embodiment for specific implementation, and will not be repeated here.

[0083] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0084] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. The apparatus embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. Furthermore, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Additionally, the shown or discussed mutual couplings, direct couplings, or communication connections may be through some communication interfaces; indirect couplings or communication connections between devices or units may be electrical, mechanical, or other forms.

[0085] In addition, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0086] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a processor-executable, non-volatile, computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0087] In this specification, the same or similar parts between the various embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, for the product embodiments described later, since they correspond to the methods, the descriptions are relatively simple, and the relevant parts can be referred to the descriptions in the system embodiments.

[0088] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for estimating motion vectors in ultrasound imaging, characterized in that, include: Acquire ultrasonic echo complex signals or their envelope sequences arranged in a time frame sequence as the initial signal; The initial signal is subjected to time high-pass filtering and spatial difference Gaussian band-pass filtering in sequence to suppress static background and enhance the spatial contrast of moving targets; Multiple Reichardt correlators in multiple directions are constructed in the IQ complex domain. Each Reichardt correlator is subjected to signal delay, conjugate multiplication, and antisymmetric differential operation to obtain multiple correlation energies in the multiple directions. A multi-scale parameter group is set, and the preferred velocity for each scale is obtained based on the pixel physical spacing and frame rate. The displacement under each scale is mapped to the preferred velocity, and the amplitude of the corresponding preferred velocity is fused with the relevant energy as weight to output an omnidirectional velocity estimate. The omnidirectional velocity estimate includes velocity amplitude and direction unit vector. The estimated axial velocity of the initial signal is obtained based on the Kasai autocorrelation method. The omnidirectional velocity estimate is adaptively weighted and fused with the axial velocity, and the motion vector and confidence level are output. The fusion weight is determined based on the Kasai autocorrelation amplitude and the correlation energy.

2. The motion vector estimation method according to claim 1, characterized in that, After performing time-pass filtering on the initial signal, an adaptive threshold is determined based on the interquartile range of the local intensity change and the median of the local intensity of the filtered initial signal. The inter-frame intensity difference of the initial signal is decomposed into a signal enhancement channel and a signal attenuation channel, and a dual-channel gating weight is constructed to suppress pseudo-intensity changes caused by coherent speckle. The adaptive threshold τ = κ·IQR_local + β·median, where κ and β are adjustment coefficients, and their values ​​range from κ∈[0.02, 0.15] and β∈[0, 0.2], respectively.

3. The motion vector estimation method according to claim 1, characterized in that, The time-pass filtering is achieved using the formula Y(t)=X(t)-αX(t-1); Where X(t) is the current frame signal, X(t-1) is the previous frame signal, and α is the filter coefficient, which takes values ​​in the range of [0.85, 0.98].

4. The motion vector estimation method according to claim 1, characterized in that, The spatial difference Gaussian bandpass filter is achieved by the difference between two Gaussian kernels with standard deviations of σ1 and σ2, where σ1 ∈ [0.4, 0.8] pixels, σ2 ∈ [1.0, 1.8] pixels, and σ1 < σ2.

5. The motion vector estimation method according to claim 4, characterized in that, The spatial difference Gaussian bandpass filter is implemented using a split convolution method, which decomposes the two-dimensional Gaussian convolution into sequentially executed horizontal and vertical one-dimensional convolutions.

6. The motion vector estimation method according to claim 1, characterized in that, The multiple directions include four directions: left, right, up, and down; among them, the right-direction related energy is calculated using the following formula: R→=Re{X(x,y,t)·X*(x+Δx,y,t+Δt)-X(x+Δx,y,t)·X*(x,y,t+Δt)}, Where X() is the signal, X*() is its conjugate, and Re represents taking the real part; the correlation energy in the left, up, and down directions is calculated by symmetrically adjusting the displacement parameters (Δx, Δy).

7. The motion vector estimation method according to claim 1, characterized in that, The multi-scale parameter set includes the number of horizontal displacement pixels Δx, the number of vertical displacement pixels Δy, and the number of time delay frames Δt; the preferred velocity v_pref is calculated using the following formula: ； Where Δ_x and Δ_y are the physical pixel spacing in the horizontal and vertical directions, respectively, and fps is the imaging frame rate.

8. The motion vector estimation method according to claim 1, characterized in that, Before adaptively weighting and fusing the omnidirectional velocity estimate and the axial velocity, the two are first aligned in terms of time length. The alignment time window length is T_align=min(T_corr, T_dop), where T_corr is the time length used to calculate the omnidirectional velocity estimate and T_dop is the time length used to calculate the axial velocity.

9. The motion vector estimation method according to claim 8, characterized in that, In the adaptive weighted fusion, the fusion weight w_d assigned to the axial velocity is determined by the following function: w_d=γ·C_d / (γ·C_d+(1-γ)·C_c); Where C_d is the normalized Kasai autocorrelation magnitude, C_c is the normalized correlation energy, and γ is an adjustable parameter with a value range of [0.3, 0.8].

10. A motion vector estimation system for ultrasound imaging, characterized in that, include: The time-domain gating module is used to acquire ultrasonic echo complex signals or their envelope sequences arranged in a time frame sequence as the initial signal; The initial signal is sequentially subjected to time-pass filtering to suppress static background; The spatial bandpass module is used to perform spatial difference Gaussian bandpass filtering on the time-domain processed signal received from the time-domain gating module to enhance the spatial contrast of the moving target. An IQ domain directional correlator array is used to construct multiple directional Reichardt correlators in the IQ complex domain. Each Reichardt correlator is processed by signal delay, conjugate multiplication, and antisymmetric differential operation to obtain multiple correlation energies in the multiple directions. The multi-scale velocity calibration and fusion module is used to receive the correlation energy output by the IQ domain direction correlator array, set a multi-scale parameter group, and obtain the preferred velocity for each scale based on the pixel physical spacing and frame rate. The displacement under each scale is mapped to the preferred velocity, and the corresponding preferred velocity amplitude is fused with the correlation energy as weight to output an omnidirectional velocity estimate. The omnidirectional velocity estimate includes velocity amplitude and direction unit vector. The Kasai Doppler module is used to receive the initial signal, calculate and output the estimated axial velocity using the Kasai autocorrelation method; The adaptive fusion and confidence module is used to receive the estimated axial velocity and perform adaptive weighted fusion of the omnidirectional velocity estimate and the axial velocity, wherein the fusion weight is determined based on the Kasai autocorrelation amplitude and the correlation energy, and finally outputs the motion vector and confidence.

11. The motion vector estimation system according to claim 10, characterized in that, The spatial bandpass module employs a split convolution method to achieve spatial difference Gaussian bandpass filtering; or The IQ domain directional correlator array consists of programmable delay line units and complex multiplier / accumulator units. The programmable delay line units perform signal delay operations, and the complex multiplier / accumulator units perform numerical calculations of conjugate multiplication and antisymmetric difference.

12. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, performs the steps of the motion vector estimation method as described in any one of claims 1 to 9.

13. An electronic device, characterized in that, include: The device includes a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the electronic device is running, the processor communicates with the memory via the bus. The machine-readable instructions are executed by the processor to perform the steps of the motion vector estimation method as described in any one of claims 1 to 9.

Images