A method, system, medium, and apparatus for motion vector estimation for ultrasound imaging

Abstract

The application discloses an ultrasonic imaging motion vector estimation method, system, medium and equipment, the method comprises sequentially performing time high-pass filtering and spatial difference Gaussian band-pass filtering on an initial signal; a plurality of directional Reichardt correlators are constructed in an IQ complex domain, and a plurality of correlation energies in multiple directions are obtained; a plurality of multi-scale parameter groups are set, and a preferred speed of each scale is obtained based on a pixel physical pitch and a frame rate; the displacement under each scale is mapped into the preferred speed, and the correlation energy is used as a weight to perform fusion, and omnidirectional velocity estimation is output, the omnidirectional velocity estimation comprises a velocity amplitude and a directional unit vector; the estimated axial velocity of the initial signal is obtained based on a Kasai autocorrelation method; the omnidirectional velocity estimation is adaptively weighted and fused with the axial velocity, and a motion vector and a confidence degree are output. The application improves the estimation performance of the motion vector, and solves the problem that the prior art is not sensitive to low-speed weak flow signals.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of ultrasonic imaging and signal processing, and in particular to an ultrasonic imaging motion vector estimation method, system, medium and device. BACKGROUND

[0002] In the field of medical ultrasonic imaging, motion vector estimation is crucial for dynamic analysis of tissue and blood flow. Currently, the mainstream clinical velocity estimation techniques mainly include the following categories: one is Doppler technology, such as color Doppler and pulse Doppler, which estimates the velocity component along the ultrasonic beam axis by the phase change of the echo signal; the other is the video processing method based on B-mode image, such as block matching or optical flow method, which infers the motion by analyzing the gray level change of pixels or regions between consecutive frames.

[0003] In the prior art, the Doppler technology is limited by its physical principle, and is very sensitive to the angle between the motion direction and the beam direction, resulting in serious underestimation of the non-axial motion component. At the same time, the pulse repetition frequency (PRF) and Nyquist limit make the upper limit of the velocity limited by PRF / Nyquist. The insensitivity to low-speed weak flow is mainly caused by factors such as wall filtering and low SNR. On the other hand, the B-mode video-based method is easily disturbed by the pseudo "flicker" generated by the inherent speckle noise of ultrasound, and false motion vectors (pseudo flow) are easily generated under low signal-to-noise ratio conditions. Although multi-angle transmission or vector flow imaging (VFI) technology can improve the angle dependence problem to some extent, these schemes usually require major modifications to the ultrasonic equipment and beam sequence, and require high system computing resources, making it difficult to achieve low-cost, high-frame-rate real-time imaging in a clinical environment.

[0004] Therefore, there is a need for a new technical solution. SUMMARY

[0005] Therefore, the embodiments of the present application provide an ultrasonic imaging motion vector estimation method, system, medium and device to at least solve the problem of insensitivity to low-speed weak flow signals in the prior art.

[0006] The embodiments of the present application provide the following technical solutions:

[0007] The embodiments of the present application provide an ultrasonic imaging motion vector estimation method, comprising:

[0008] Acquire an ultrasonic echo complex signal or an envelope sequence arranged in a time frame sequence as an initial signal;

[0009] 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 the motion target;

[0010] constructing multiple directional Reichardt correlators in IQ complex domain, each of the Reichardt correlators being operated by signal delay, conjugate multiplication and antisymmetric difference, and obtaining multiple correlation energies of the multiple directions;

[0011] setting a multi-scale parameter group, and obtaining a preferred velocity of each scale based on pixel physical pitch and frame rate, mapping displacement at each scale to the preferred velocity, fusing amplitudes of the preferred velocities corresponding to the correlation energies as weights, and outputting an omnidirectional velocity estimation including a velocity amplitude and a directional unit vector;

[0012] obtaining an estimated axial velocity of the initial signal based on Kasai autocorrelation method, adaptively weighting and fusing the omnidirectional velocity estimation and the axial velocity, and outputting a motion vector and a confidence level, wherein a weighting weight is determined based on a Kasai autocorrelation amplitude and the correlation energy.

[0013] Further, after the time high-pass filtering on the initial signal, an adaptive threshold is determined based on a local intensity variation quartile range and a local intensity median of the filtered initial signal, an inter-frame intensity difference of the initial signal is decomposed into a signal enhancement channel and a signal weakening channel, and a dual-channel gating weight is constructed to suppress false intensity variation caused by coherent speckles, wherein the adaptive threshold is τ = κ · IQR_local + β · median, and κ and β are adjustment coefficients, and the value ranges of κ and β are κ ∈ [0.02, 0.15] and β ∈ [0, 0.2] respectively.

[0014] Further, the time high-pass filtering is implemented by the formula Y(t) = X(t) - αX(t-1);

[0015] wherein X(t) is a current frame signal, X(t-1) is a previous frame signal, and α is a filtering coefficient, and the value range of α is [0.85, 0.98].

[0016] Further, the spatial difference Gaussian band-pass filtering is implemented by a difference value of two Gaussian kernels with standard deviations σ1 and σ2, wherein σ1 ∈ [0.4, 0.8] pixels, σ2 ∈ [1.0, 1.8] pixels, and σ1 < σ2.

[0017] Further, the spatial difference Gaussian band-pass filtering is implemented in a separate convolution manner, and two-dimensional Gaussian convolution is decomposed into sequentially executed horizontal and vertical one-dimensional convolution.

[0018] Further, the multiple directions include left, right, up and down directions; wherein the right direction correlation energy is calculated by the following formula:

[0019] R→=Re{X(x, y, t)·X*(x+Δx, y, t+Δt)-X(x+Δx, y, t)·X*(x, y, t+Δt)},

[0020] where X() is the signal, X*() is its conjugate, and Re represents taking the real part; the correlation energy in the left, upper, and lower directions is calculated by symmetrically adjusting the displacement parameters (Δx, Δy).

[0021] Further, the multi-scale parameter set includes the horizontal displacement pixel number Δx, the vertical displacement pixel number Δy, and the time delay frame number Δt; the preferred speed v_pref is calculated by the following formula:

[0022] v_pref=√((Δx·Δ_x)^2+(Δy·Δ_y)^2) / (Δt / fps);

[0023] where Δ_x and Δ_y are the physical pixel pitches in the horizontal and vertical directions, respectively, and fps is the imaging frame rate.

[0024] Further, before adaptively and weightedly fusing the omnidirectional velocity estimation and the axial velocity, the two are first time length aligned, and the aligned time window length T_align=min(T_corr, T_dop), where T_corr is the time length used to calculate the omnidirectional velocity estimation, and T_dop is the time length used to calculate the axial velocity.

[0025] Further, in the adaptive weighted fusion, the fusion weight w_d given to the axial velocity is determined by the following function:

[0026] w_d=γ·C_d / (γ·C_d+(1-γ)·C_c);

[0027] where C_d is the normalized Kasai autocorrelation amplitude, C_c is the normalized correlation energy, and γ is an adjustable parameter, with a value range of [0.3, 0.8].

[0028] The present application provides an ultrasonic imaging motion vector estimation system, comprising:

[0029] a time domain gating module, configured to collect ultrasonic echo complex signals or envelope sequences arranged in time frame sequences as initial signals; and sequentially perform time high-pass filtering on the initial signals to suppress static background;

[0030] a spatial bandpass module, configured to perform spatial difference Gaussian bandpass filtering on the time domain processed signals received from the time domain gating module to enhance the spatial contrast of the moving target;

[0031] An IQ domain direction correlator array is configured to construct a plurality of direction Reichardt correlators in an IQ complex domain, each of the Reichardt correlators is obtained by signal delay, conjugate multiplication and anti-symmetry difference operation, and a plurality of correlation energies of the plurality of directions are obtained;

[0032] A multi-scale velocity calibration and fusion module is configured to receive the correlation energies output by the IQ domain direction correlator array, set a multi-scale parameter group, obtain a preferred velocity of each scale based on a pixel physical pitch and a frame rate, map a displacement under each scale to the preferred velocity, fuse magnitudes of the preferred velocities corresponding to the correlation energies as weights, and output an omnidirectional velocity estimation, the omnidirectional velocity estimation including a velocity magnitude and a direction unit vector;

[0033] A Kasai Doppler module is configured to receive the initial signal, calculate and output an estimated axial velocity by using a Kasai autocorrelation method;

[0034] An adaptive fusion and confidence module is configured to receive the estimated axial velocity, and adaptively and weightedly fuse the omnidirectional velocity estimation and the axial velocity, wherein a fusion weight is determined based on a Kasai autocorrelation magnitude and the correlation energy, and finally output a motion vector and a confidence.

[0035] Further, the spatial band-pass module is configured to realize spatial difference Gaussian band-pass filtering in a separate convolution manner; or

[0036] The IQ domain direction correlator array is composed of a programmable delay line unit and a complex multiplication / accumulation unit, the programmable delay line unit is configured to perform signal delay operation, and the complex multiplication / accumulation unit is configured to perform numerical calculation of conjugate multiplication and anti-symmetry difference.

[0037] The application further provides a computer readable storage medium, the computer readable storage medium stores a computer program, and the computer program is executed by a processor to perform the steps of the motion vector estimation method.

[0038] The application further provides an electronic device, which comprises a processor, a memory and a bus, the memory stores machine readable instructions executable by the processor, the processor and the memory communicate through the bus when the electronic device is running, and the machine readable instructions are executed by the processor to perform the steps of the motion vector estimation method.

[0039] Compared with the prior art, the above at least one technical solution adopted by the embodiment of the application can achieve at least the following beneficial effects:

[0040] The application improves the estimation performance of the motion vector through the synergistic technical chain of the space-time pretreatment-all direction velocity extraction-axis velocity fusion, improves the sensitivity to weak flow and slow flow through the time high-pass and space DoG band-pass filtering, further improves the non-axis detection capability through the IQ complex domain four-direction Reichardt correlator, and uses the complex signal phase information to completely capture the transverse motion component. BRIEF DESCRIPTION OF DRAWINGS

[0041] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed to be used in the embodiments will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0042] Fig. 1-2 A flow chart of an ultrasound imaging motion vector estimation method according to an embodiment of the present application;

[0043] Fig. 3 A time high-pass filtering and IQR-based ON / OFF gating timing diagram according to an embodiment of the present application;

[0044] Fig. 4 An IQ domain four-direction Reichardt correlator structure and signal path diagram according to an embodiment of the present application;

[0045] Fig. 5 A Kasai autocorrelation and adaptive fusion weight function diagram of correlator energy according to an embodiment of the present application;

[0046] Fig. 6 An ON / OFF gating flow chart according to an embodiment of the present application;

[0047] Fig. 7 A velocity / direction superimposed diagram and confidence diagram of a software embodiment according to an embodiment of the present application;

[0048] Fig. 8 ROC / AUC and minimum detectable speed (MDS) curves of a comparison and ablation experiment according to an embodiment of the present application. DETAILED DESCRIPTION

[0049] The embodiments of the present application will be described in detail below with reference to the drawings.

[0050] 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.

[0051] 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.

[0052] 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.

[0053] 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.

[0054] 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.

[0055] Based on this, the embodiments of this specification propose a processing solution: such as Fig. 1 As 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.

[0056] Specifically, the motion vector estimation method of the present invention can be summarized as a cat's-eye spatiotemporal signal chain:

[0057] Time-based high-pass filtering + ON / OFF gating: First-order difference high-pass filtering of the frame sequence Y(t) = X(t) αX(t 1), and based on the adaptive threshold of local IQR, the intensity change is decomposed into ON (enhancement) and OFF (decrease) channels to form gating weights to suppress speckle spurious flicker;

[0058] Spatial bandpass (DoG): Enhances the spatial contrast between moving edges and speckle clusters;

[0059] 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;

[0060] 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;

[0061] 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;

[0062] Time length alignment and post-processing: T_align=min(T_corr,T_dop), outputs velocity vector, direction angle and confidence level.

[0063] The technical scheme provided by the embodiments of the present application is described below with reference to the drawings.

[0064] As shown in Fig. 1-2 The present application provides a motion vector estimation method for ultrasonic imaging, comprising:

[0065] In step S102, an ultrasonic echo complex signal or an envelope sequence thereof arranged in a time frame sequence is collected as an initial signal.

[0066] The ultrasonic echo complex signal (IQ signal) is a complex signal composed of a real part (I) and an imaginary part (Q). It can be expressed as A·e^(jφ), where A is the amplitude and φ is the phase. In an ultrasonic system, it is obtained after quadrature demodulation of a 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 number form, i.e. the B-mode (brightness mode) grayscale image we usually see.

[0067] Wherein, the ultrasonic probe can be used to emit ultrasonic waves into the human body and receive echo signals scattered or reflected back from tissues and blood flow.

[0068] Wherein, the initial signals are arranged one by one in time sequence to form a three-dimensional data set.

[0069] Step S102 is used to provide the original input required by all subsequent processing modules (time domain filtering, spatial filtering, correlator, Doppler estimation).

[0070] In step S104, the initial signal is sequentially subjected to time high-pass filtering and spatial difference Gaussian band-pass filtering to suppress static background and enhance the spatial contrast of the moving target.

[0071] As shown in Fig. 3 In the step, the initial signal is first subjected to time high-pass filtering, and then the pseudo signal is processed and subjected to spatial difference Gaussian band-pass filtering.

[0072] Wherein, the time high-pass filtering is realized by the formula Y(t)=X(t)-αX(t-1); wherein, X(t) is the current frame signal, X(t-1) is the previous frame signal, and α is the filtering coefficient, which is in the range of [0.85, 0.98].

[0073] Specifically, since the static or slowly changing background tissue is very similar in adjacent frames, the result after subtraction is close to zero, so that the frames of the static background are filtered out by the above formula; since the moving target changes in position or form between frames, a significant difference signal is generated after subtraction, so that it is highlighted.

[0074] Therefore, the temporal high-pass filtering is used to filter out the static background signal which does not change over time, and provide a "clean" time series signal containing only motion information for the subsequent motion estimation step.

[0075] After the initial signal is temporally high-pass filtered, the method further comprises the steps of:

[0076] In step S104a, based on the local intensity change interquartile range and the local intensity median 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 weakening channel, and a double-channel gating weight is constructed to suppress the false intensity change caused by coherent speckle, wherein the adaptive threshold is as follows:

[0077] τ=κ·IQR_local+β·median, κ∈[0.02, 0.15], β∈[0, 0.2];

[0078] Wherein, κ and β are adjustment coefficients, and the value ranges thereof are κ ∈ [0.02, 0.15] and β ∈ [0, 0.2] respectively.

[0079] Wherein, the local intensity change interquartile range is used to measure the dispersion degree of signal change in the local area; the local intensity median represents the reference level of the signal intensity in the area; the adaptive threshold is used to provide a dynamic judgment standard related to the statistical characteristics of the local image, avoiding the problem that the fixed threshold is too sensitive (introducing noise) in some areas or not sensitive (losing real signal) in other areas.

[0080] Wherein, the signal enhancement channel only retains the part where the inter-frame intensity difference is greater than zero, and the signal weakening channel only retains the part where the inter-frame intensity difference is less than zero, and the comparison result of 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, it will be given a higher weight.

[0081] Wherein, the gating weight is used to regulate the "trust degree" of the current pixel signal in the subsequent processing step. The weight of 1 means that the change at this place is completely trusted to be real motion; the weight of 0 means that the change at this place is completely ignored and considered as noise; and the weight between them means partial trust. This weight is directly multiplied by the signal of the point, thereby suppressing noise at the source.

[0082] Wherein, step S104a is located after the temporal high-pass filtering, and is used to distinguish the real motion signal and the false intensity change caused by coherent speckle noise, and by processing the positive and negative signal changes through the double gates respectively, the system can process complex motion patterns more finely.

[0083] The spatial difference Gaussian band-pass filtering is realized 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.

[0084] The spatial difference Gaussian band-pass filtering is realized 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.

[0085] Specifically, the small standard deviation Gaussian kernel (σ1) is used to preserve small edges and details (possibly containing noise), while the large standard deviation Gaussian kernel (σ2) is used to capture large-scale background structures. Subtracting the two, the large-scale uniform background is offset, while the features between the two scales are preserved and enhanced.

[0086] The spatial difference Gaussian band-pass filtering weakens the uniform region in the spatial dimension and sharpens the boundary and internal texture structure of the moving target, thereby improving the spatial distinguishability of the moving target.

[0087] Further, the spatial difference Gaussian band-pass filtering is realized in a separate convolution manner, which decomposes the two-dimensional Gaussian convolution into sequentially executed horizontal and vertical one-dimensional convolutions.

[0088] The step S104 realizes the "time domain first, then spatial domain" cooperative preprocessing by sequentially performing the time high-pass filtering and the 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 features of the moving target, so that the contrast between the moving target and its background in the output signal is greatly improved. The IQ domain direction correlator is provided with high-quality and high signal-to-noise ratio input, which directly improves the detection ability and estimation accuracy of the entire system for weak and non-axial motion.

[0089] The step S106, in the IQ complex domain, constructs a plurality of direction Reichardt correlators, each of which is obtained by signal delay, conjugate multiplication and anti-symmetric difference operation, and obtains a plurality of correlation energies in multiple directions.

[0090] In combination with Fig. 4As shown, the Reichardt correlator is a motion detection model that mainly perceives motion in a specific direction by comparing the signals of two adjacent spatial points at different times. The signal delay is a small time delay At to the signal of one of the spatial points, which simulates the time required for the motion signal to propagate from one point to another; the conjugate multiplication is to multiply the undelayed signal and the conjugate of the delayed signal by a complex number, which is equivalent to calculating the cross-correlation of the two signals in signal processing, and the result is a complex number, whose phase contains the phase difference information between the two signals, and the phase difference is extremely sensitive to motion; the anti-symmetry difference calculates the difference between the correlation results of two opposite directions (for example, right correlation minus left correlation). This is to give the response direction selectivity.

[0091] wherein the plurality of directions include left, right, up, and down directions; wherein the right correlation energy is calculated by the following formula:

[0092] R→=Re{X(x,y,t)·X*(x+Δx,y,t+Δt)-X(x+Δx,y,t)·X*(x,y,t+Δt)},

[0093] wherein X() is a signal, X*() is its conjugate, and Re represents taking the real part; the correlation energies of the left, up, and down directions are calculated by symmetrically adjusting the displacement parameters (Δx, Δy).

[0094] Step S106 is used to detect and quantify the motion intensity (correlation energy) in the four basic directions of up, down, left, and right. The Reichardt correlator in the four directions can perceive and measure the transverse motion component perpendicular to the direction of the ultrasonic beam, thereby overcoming the fundamental limitation of the conventional Doppler technique that can only measure the axial velocity. In addition, step S106 outputs a clear "correlation energy" value for each direction. The size of this value represents the intensity of the motion signal in that direction, and the positive or negative value represents the direction of the motion. This lays the foundation for generating an accurate vector direction in subsequent fusion.

[0095] Step S106 can utilize phase information by operating in the IQ complex domain, making it very sensitive to extremely slow motion (slow flow) and real motion under low signal-to-noise ratio conditions, and the effect is much better than the traditional optical flow method running on B-mode images.

[0096] Step S108, set a plurality of scale parameter groups, and based on the physical spacing of the pixels and the frame rate, map the displacement under each scale to the preferred velocity, and fuse the corresponding preferred velocity amplitudes with the correlation energy as the weight, and output the omnidirectional velocity estimation, which includes the velocity amplitude and the direction unit vector.

[0097] wherein the multi-scale parameter set includes horizontal displacement pixel number Δx, vertical displacement pixel number Δy, and time delay frame number Δt; the preferred velocity v pref is calculated by the following formula:

[0098] v pref = √((Δx·Δ x)2+(Δy·Δ y)2) / (Δt / fps);

[0099] wherein Δ x and Δ y are respectively the horizontal and vertical pixel physical pitch, and fps is the imaging frame rate.

[0100] wherein the multi-scale parameter set (Δx, Δy, Δt) is a series of preset combinations covering different displacement ranges and time delays. For example, one set of parameters can be for large displacement fast flow, and another set for small displacement slow flow.

[0101] Specifically, a multi-scale parameter set (including horizontal displacement pixel number Δx, vertical displacement pixel number Δy, and time delay frame number Δt) is set, and based on the pixel physical pitch (horizontal Δ x, vertical Δ y) and the frame rate (fps), the displacement under each scale is mapped to the preferred velocity of the corresponding scale: wherein the preferred velocity amplitude of each scale is calculated by the formula v pref = √((Δx·Δ x)2+(Δy·Δ y)2) / (Δt / fps), and the direction vector is determined by the normalized displacement direction (Δx, Δy) of the scale; then the preferred velocity amplitudes of all scales are weighted and summed (i.e. the final velocity amplitude = Σ (correlation energy i × v pref i) / Σ correlation energy i) with the relevant energy of each scale as the weight, and the direction unit vectors of all scales are weighted and summed and then normalized again (i.e. the final direction unit vector = [Σ (correlation energy i × direction vector i)] / Σ correlation energy i), and the omnidirectional velocity estimation is output by the above fusion, which includes the velocity amplitude and the direction unit vector.

[0102] wherein the omnidirectional velocity estimation is the final output of this step, which integrates information from all directions and all scales to represent the motion state of each pixel point in a single vector form. The vector is composed of two parts: the velocity amplitude (scalar, indicating speed) and the direction unit vector (only indicating direction, length 1), which completely describes the speed and direction of motion.

[0103] Step S106 detects by multi-scale, which can ensure that the system can simultaneously sensitively detect fast flow and slow flow, and through energy weighted fusion, it can also automatically filter and fuse the most reliable and consistent final estimate value from multiple tentative velocities, thereby improving the accuracy and robustness of the estimate.

[0104] In step S110, the estimated axial velocity of the initial signal is obtained based on the Kasai autocorrelation method, the omnidirectional velocity estimation and the axial velocity are adaptively weighted and fused, and the motion vector and the confidence are output, wherein the fusion weight is determined based on the Kasai autocorrelation amplitude and the correlation energy.

[0105] The Kasai autocorrelation method is a Doppler velocity estimation algorithm, which estimates the phase change rate by calculating the autocorrelation function between adjacent ultrasonic echo signals, so as to calculate the velocity component along the ultrasonic beam direction. Its advantages are high calculation efficiency, accurate axial velocity measurement, and is the standard algorithm for color Doppler imaging.

[0106] Before adaptively weighting and fusing the omnidirectional velocity estimation and the axial velocity, the time length of the two is aligned, and the aligned time window length T align = min (T corr, T dop), wherein T corr is the time length used for calculating the omnidirectional velocity estimation, and T dop is the time length used for calculating the axial velocity.

[0107] The time window length formula determines the common effective time window, and realizes synchronization by truncating the longer data sequence, thereby preventing errors caused by time misalignment.

[0108] Further, in the adaptive weighted fusion, the fusion weight w d of the axial velocity is determined by the following function:

[0109] w d = γ·C d / (γ·C d + (1-γ)·C c) ;

[0110] Wherein, C d is the normalized Kasai autocorrelation amplitude, which is the confidence of the axial velocity path, C c is the normalized correlation energy, which is the confidence of the omnidirectional velocity path, and γ is an adjustable parameter, which is in the range of [0.3, 0.8], and is used to identify which path is biased as a whole.

[0111] As shown in FIG. Fig. 5-8 Through time alignment, step S110 can ensure the time consistency of fusion, which is used to ensure the smoothness and accuracy of the final motion vector in time, and the adaptive weight is used to dynamically determine which path the final result depends on, so that the optimal estimation can be obtained in different regions such as the center of the blood vessel (axial flow dominant) and the edge / curve (lateral flow significant), thereby outputting the optimal motion vector.

[0112] The motion vector estimation method has the following beneficial effects: weak flow / slow flow sensitive, strong non-axial ability, high robustness (IQ domain correlation / conjugate multiplication), resource friendly (DoG separation convolution, delay line+complex multiplication MAC pipeline), and strong system-level availability (can be deployed in the form of software / firmware upgrade on the existing system).

[0113] The motion vector estimation system for ultrasonic imaging has the axial precision of the Kasai method and the lateral sensitivity of the Reichardt correlator, and can output a confidence map quantifying the reliability of the velocity estimation value of each pixel point, so that doctors can determine which blood flow signals are reliable and which may have uncertainty, and the reliability of clinical diagnosis is greatly improved.

[0114] The application further provides a motion vector estimation system for ultrasonic imaging, which comprises a time domain gating module, a spatial band-pass 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 data with the multi-scale velocity calibration and fusion module and an adaptive fusion and confidence module.

[0115] The time domain gating module is used for collecting ultrasonic echo complex signals or envelope sequences arranged in time frame sequences as initial signals, performing time high-pass filtering on the initial signals in sequence to suppress static background, determining an adaptive threshold based on the local intensity variation quartile range and local intensity median of the filtered signals, decomposing the inter-frame intensity difference into a signal enhancement channel and a signal weakening channel, constructing a double-channel gating weight to suppress the false intensity variation caused by coherent speckle, and outputting a time domain processed signal to the spatial band-pass module.

[0116] The spatial band-pass module is used for performing spatial difference Gaussian band-pass filtering on the time domain processed signal output by the time domain gating module to enhance the spatial contrast of the moving target.

[0117] Specifically, the spatial band-pass module performs spatial difference Gaussian band-pass filtering through the difference between two Gaussian kernels with standard deviations σ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 a spatial processed signal to the IQ domain direction correlator array.

[0118] The IQ domain direction correlator array is used for constructing a plurality of direction Reichardt correlators in the IQ complex domain, each Reichardt correlator being obtained through signal delay, conjugate multiplication and antisymmetric difference operation, and obtaining a plurality of correlation energies in a plurality of directions.

[0119] Specifically, the IQ domain direction correlator array receives the spatial processed signals output by the spatial band-pass module, constructs left, right, top and bottom four directions of Reichardt correlators in the IQ complex domain, each of the Reichardt correlators sequentially performs signal delay, conjugate multiplication and anti-symmetry difference operation, calculates and outputs the correlation energy of the four directions, and sends the correlation energy to the multi-scale velocity calibration and fusion module and the adaptive fusion and confidence module respectively.

[0120] The multi-scale velocity calibration and fusion module is configured to receive the correlation energy output by the IQ domain direction correlator array, set a multi-scale parameter group, and obtain a preferred velocity of each scale based on a pixel physical pitch and a frame rate, map a displacement under each scale to the preferred velocity, fuse the preferred velocities under each scale by using the correlation energy as a weight, and output an omnidirectional velocity estimation including a velocity amplitude and a direction unit vector.

[0121] 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 group including a horizontal displacement pixel number Δx, a vertical displacement pixel number Δy and a time delay frame number Δt, calculates a preferred velocity of each scale according to a preferred velocity formula based on a pixel physical pitch (Δ_x, Δ_y) and an imaging frame rate (fps) in the horizontal and vertical directions, fuses the preferred velocities of each scale by using the correlation energy as a weight, and outputs the omnidirectional velocity estimation including the velocity amplitude and the direction unit vector to the adaptive fusion and confidence module.

[0122] The Kasai Doppler module is configured to receive the initial signal and calculate and output an estimated axial velocity by using a Kasai autocorrelation method.

[0123] The adaptive fusion and confidence module is configured to receive the estimated axial velocity, and adaptively and weightedly fuse the omnidirectional velocity estimation and the axial velocity, wherein a fusion weight is determined based on a Kasai autocorrelation amplitude and the correlation energy, and finally outputs a motion vector and a confidence.

[0124] Specifically, the adaptive fusion and confidence module receives the omnidirectional velocity estimation output by the multi-scale velocity calibration and fusion module, a time length T_corr of omnidirectional velocity calculation, the estimated axial velocity output by the Kasai Doppler module, a time length T_dop of axial velocity calculation, and the correlation energy output by the IQ domain direction correlator array, first performs time length alignment on the omnidirectional velocity estimation and the estimated axial velocity according to T_align = min (T_corr, T_dop), and then based on a Kasai autocorrelation amplitude (normalized as C_d) and the correlation energy (normalized as C_c), performs fusion according to a formula

[0125] The fusion weight is determined, and adaptive weighted fusion is performed on the aligned omnidirectional velocity estimation and the estimated axial velocity, and finally the motion vector and the confidence are output.

[0126] Further, the spatial band-pass module implements the spatial difference Gaussian band-pass filtering in a separated convolution manner, decomposes the two-dimensional Gaussian convolution into sequentially executed horizontal one-dimensional convolution and vertical one-dimensional convolution, and reduces the calculation complexity to adapt to real-time processing.

[0127] Further, the IQ domain direction correlator array is composed of a programmable delay line unit and a complex multiplication / accumulation unit, the programmable delay line unit performs a signal delay operation, and the complex multiplication / accumulation unit performs conjugate multiplication and numerical calculation of anti-symmetrical difference, adapting to the pipeline data processing of FPGA / SoC.

[0128] The embodiment of the application further provides a computer readable storage medium, which stores a computer program, and the computer program can execute the steps of the motion vector estimation method in the method embodiment when the computer program is run by a processor. Fig. 1 The specific implementation can refer to the method embodiment, and will not be described here.

[0129] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the system, device and unit described above can refer to the corresponding process in the foregoing method embodiment, and will not be described here.

[0130] In several embodiments provided in the application, it should be understood that the disclosed system, device and method can be implemented in other ways. The device embodiments described above are only schematic, for example, the division of the units is only a logical function division, and actual implementation can have another division manner, for example, a plurality of units or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the units shown or discussed can be indirect coupling or communication connection through some communication interface, device or unit, and can be electrical, mechanical or other forms.

[0131] In addition, each functional unit in each embodiment of the application can be integrated in one processing unit, or each unit can exist physically, or two or more units can be integrated in one unit.

[0132] If the functions are implemented in the form of software function units and sold or used as independent products, they can be stored in a nonvolatile computer readable storage medium executable by a processor. Based on this understanding, the technical solutions of the present application essentially or the parts that contribute to the prior art or parts of the technical solutions can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes a number of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present application. The aforementioned storage medium includes: a U disk, a mobile hard disk, a read-only memory (Read-Only Memory, ROM), a random access memory (Random Access Memory, RAM), a magnetic disk or an optical disk, and various media that can store program codes.

[0133] In the specification, the same or similar parts between various embodiments can be referred to each other, and each embodiment focuses on the difference from other embodiments. Especially, for the product embodiment described later, since it is corresponding to the method, the description is relatively simple, and the relevant part can be referred to the part of the system embodiment.

[0134] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited to this. Any person skilled in the art can easily think of changes or replacements within the technical scope disclosed in the present application, which should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection 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 the time-pass high-pass filtering on the initial signal, an adaptive threshold is determined based on the interquartile range of the local intensity variation and the median 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 variations caused by coherent speckle. The adaptive threshold τ = κ·IQR_local + β·median, where κ and β are adjustment coefficients with values ​​ranging from κ∈[0.02, 0.15] and β∈[0, 0.2], respectively. IQR_local is the interquartile range of the local intensity variation, and median is the median local intensity.

3. The motion vector estimation method according to claim 1, characterized in that, 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 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: v_pref=√((Δx·Δ_x)^2+(Δy·Δ_y)^2) / (Δt / fps); 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