CN114240815A - Multithreading strain imaging method and device based on living body ultrasonic image - Google Patents

Abstract

The invention discloses a multithreading strain imaging method and a multithreading strain imaging device based on a living body ultrasonic image, and belongs to the field of medical image processing. Collecting a living body B ultrasonic image sequence, analyzing the correlation characteristics of the images, judging the respiratory frequency, and dividing the images into different respiratory cycles and synchronous phase periods among the cycles; establishing a computing thread for each phase period, computing strain in parallel by each thread after displacement correction and image registration, and denoising according to the correlation among the threads; and finally, fusing the multithreading strain results to obtain accurate and high-time-density strain output. The method can solve the problems of artifacts and errors generated in strain imaging by physiological motion in the prior art, can obtain strain distribution in multiple motion phase periods, and further obtains a precise, high-robustness and high-time-density in-vivo strain distribution image after denoising and fusion are carried out by combining with the spatial correlation of the strain distribution.

Description

Multithreading strain imaging method and device based on living body ultrasonic image

Technical Field

The invention relates to the field of medical image processing, in particular to a multithreading strain imaging method and a multithreading strain imaging device based on a living body ultrasonic image.

Background

Ultrasonic imaging is one of four major image diagnosis technologies in clinical medicine, and has high application value in clinic due to the advantages of safety, no wound, small burden on patients, good signal processing real-time performance and the like. The development of ultrasonic imaging technology goes from one-dimensional to three-dimensional, static to dynamic, structural imaging to functional imaging, the function of ultrasonic diagnostic equipment is gradually enhanced, and the provided information is increasingly rich.

With the development of ultrasound imaging technology and ultrasound imaging equipment, a variety of diagnostic, analytical imaging algorithms based on ultrasound imaging signals are generated in succession. Such as doppler flow imaging, elastography, thermal strain imaging, etc., based on B-mode ultrasound imaging signals. The derived imaging algorithm quantifies the information of blood flow, hardness, temperature and the like of target tissues, greatly enhances the monitoring capability of B-ultrasonic imaging and enriches the clinical diagnosis basis.

The accuracy and signal-to-noise ratio of strain imaging based on B-mode ultrasonic imaging signals are heavily dependent on the correlation of images before and after compression. Under motion of a living body, such as breathing, heartbeat, pulse, etc., the motion causes severe artifacts in strain imaging and can even overwhelm the effective strain. Therefore, overcoming tissue motion and realizing image registration are an urgent problem to be solved in the field of strain imaging.

In the prior art, some technical solutions are proposed for imaging to overcome physiological motion of a living body. For example, the invention and creation names are: an ultrasound image generation system for generating intravascular ultrasound images (application date: 2018, 11/13; application number: 201880087709.6) discloses a method for intravascular ultrasound imaging based on intravascular pressure in the same moving state of a vessel wall. The approach uses intravascular pressure as a selection criterion for combining intravascular ultrasound images, combining image signals for which the pressure values differ from the reference pressure value by less than a predefined deviation threshold, thereby reducing artifacts due to pulse in the combined image. The disadvantage of this solution is that it requires additional blood pressure measuring equipment, and the synchronization of the signals of the manometer and the ultrasound imaging equipment and the noise signal of the manometer may cause gating failures. For another example, the invention is named as: an ultrasonic method for measuring temperature change of biological tissues based on a thermal expansion and gating algorithm (application date: 2017, 9, 25; application number: 201710876349) discloses an ultrasonic method for measuring temperature change of tissues based on a thermal expansion and gating algorithm. And selecting a matched image based on the B ultrasonic time sequence image, and training an adaptive filter by adopting a noise signal so as to remove the noise of the motion. The scheme has the defects that a mechanism for judging the physiological motion state and the motion cycle is lacked, so that accumulated motion artifacts caused by improper image matching intervals exist; and the image matching can be only carried out in the single motion state of respiration, the noise influence of the single motion state is easy to be caused, and the time density of the output strain result is low, so that the real-time requirement of clinical monitoring is not favorably met. The invention also has the name: a respiration separation type strain imaging method based on a living body ultrasound image (application date: 7/27/2020, application number: 202010731441.4) is an invention patent applied by the applicant in the earlier stage, the invention solves the problems of judging the respiration period and respectively carrying out image registration in the expiration stage and the inspiration stage, but the time interval of the output result is not uniform enough, image registration and displacement correction can not be respectively carried out in different phase stages, the space correlation between the results in different states is not utilized for carrying out strain result denoising, and a certain distance is kept away from the clinical application requirement.

In summary, it is an urgent need in the art to output strain imaging with high accuracy and meeting the clinical real-time requirement without the assistance of additional equipment, and effectively suppress the influence of motion on the image.

Disclosure of Invention

1. Problems to be solved

Aiming at the defects that motion suppression is realized through auxiliary gating of additional equipment in the prior art, and the gating technology depends on a single motion phase, the output result is low in time density and poor in real-time performance and is often accompanied with artifacts and large errors, the invention provides a multithreading strain imaging method and device based on a living body ultrasonic image, strain distribution in a multi-motion phase period can be obtained, and a living body strain distribution image with accuracy, high robustness and high time density is obtained after denoising and fusion are further combined with the spatial correlation of the strain distribution.

2. Technical scheme

In order to overcome the problems, the technical scheme provided by the invention is as follows:

a multithreading strain imaging method based on a living body ultrasonic image comprises the following steps:

collecting an ultrasonic image of a living object to obtain a digital ultrasonic image sequence I;

a characteristic extraction step: extracting the texture features of the digital ultrasonic image sequence I;

a cycle dividing step: determining characteristic images of the expiration and inspiration states through the texture characteristics of the digital ultrasonic image sequence I, and dividing the digital ultrasonic image sequence I into a plurality of image sequences of expiration periods and a plurality of image sequences of inspiration periods;

setting a threshold value: setting a expiration textural feature threshold, and if the textural feature of the image in the expiration period is greater than the expiration textural feature threshold, placing the image in an expiration image sequence I_exhDividing the image sequence of each expiratory cycle into U phases₁A phase period, U₁For the number of divided phase periods, a multithreading subset I is obtained_exh,u1(u₁＝1,2,…,U₁) (ii) a Setting an inspiration textural feature threshold, and if the textural feature of the image in the inspiration period is larger than the inspiration textural feature threshold, placing the image in an inspiration image sequence I_inhDividing the image sequence of each inspiration cycle into U phases₂A phase period, U₂For the number of divided phase periods, a multithreading subset I is obtained_inh,u2(u₂＝1,2,…,U₂)；

A registration step: determining the region of interest by taking the position of the target region as the center, for I_exh,u1、I_inh,u2To the region of interest, displacement calibration and image registration are performed to obtain registered subsets Q_exh,u1、Q_inh,u2；Q_exh,u1And Q_inh,u2Is denoted as Q_u(U-1, 2, …, U) and U is Q_uNumber of threads of (U) ═ U₁+U₂；

A fusion step: parallel computation Q_uThe time domain interpolation is carried out on the accumulated strain image sequence of each thread, denoising is carried out based on strain correlation, and the accumulated strain of each thread after processing is fused to obtain a final strain imaging result.

Furthermore, the texture features in the feature extraction step are correlations between images, and the specific extraction method comprises the following steps:

taking N-1, 2, …, M frames of N frames of image of digital ultrasonic image sequence I as template frame I_refN is the total frame number of the image sequence I, M is a positive integer, and the M frame images at least cover a complete respiration period, and the rest N-1 frames are used as detection frames I_tarCalculating I_refAnd I_tarObtaining M cross correlation coefficients-time curve p_n(t)。

Furthermore, the specific implementation method of the cycle dividing step is as follows:

within the interval of possible physiological movement frequencies, from p_n(t) frequency spectrum P_n(f) Characteristic image I for determining expiration and inspiration states_RE、I_RIAnd the respiratory frequency f_resThe calculation method comprises the following steps: will spectrum P_n(f) Is denoted as A_n，A_nThe template frame corresponding to the maximum value of (n-1, 2, …, M) is I_RECorresponding to p_n(t) is denoted by p_exh(t) for p_exh(t) after low-pass filtering, the detection frame corresponding to the minimum value in the first M points on the curve is I_RI，I_RIAs template frame correspondences_n(t) is denoted by p_inh(t)，p_exh(t) the frequency corresponding to the maximum of the spectral amplitude is f_res(ii) a To I_RECorresponding correlation coefficient-time curve p_exh(t) filtering to obtain a smooth curve p'_exhAnd (t) taking the corresponding time of all the peaks as the boundary of the inspiration period and the corresponding time of the valleys as the boundary of the expiration period, thereby dividing the digital ultrasonic image sequence I into a plurality of image sequences of expiration periods and a plurality of image sequences of inspiration periods.

Further, p 'is selected'_exhThe method of the peak and the valley of (t) is: the exhalation characteristic image I_RESet as an anchor point at a breathing frequency f_resDesign search step Δ t, 1/f_res>Δt>0.5/f_resIn the direction of increasing and decreasing time from anchor point, at p'_exh(t) alternately searching for valleys and peaks.

Furthermore, the dividing of the motion phase for the image sequence of each expiratory cycle in the setting of the threshold step is performed by dividing or setting the division of the amplitude interval of the correlation coefficient at equal time intervals.

Furthermore, the specific implementation method of the step of setting the threshold value is as follows:

at p_exh(t) in each expiratory cycle, setting the upper quartile of all cross-correlation coefficients as the expiratory texture feature threshold, for p_exh(t) points above the expiratory texture feature threshold, and placing the corresponding images into the expiratory image sequence I_exh；

Dividing the image sequence of each expiration period into U phases₁A phase period, U₁For the number of divided phase periods, a multithreading subset I is obtained_exh,u1(u₁＝1,2,…,U₁)；

At p_inh(t) in each inspiration period, setting the upper quartile of all cross-correlation coefficients as the inspiration texture feature threshold, for p_inh(t) points above the threshold of the inspiratory texture feature, and placing the corresponding images in the inspiratory image sequence I_inh；

The image sequence of each inspiration cycle is operatedDivision of the dynamic phase into U₂A phase period, U₂For the number of divided phase periods, a multithreading subset I is obtained_inh,u2(u₂＝1,2,…,U₂)。

Further, the specific implementation manner of the registration step is as follows:

the position of the target area is taken as the center, the size is determined to be AxL, and the coordinate of the upper left corner is determined to be (a)₀,l₀) Wherein a is the axial length and L is the transverse length;

during each phase (u) of the first expiratory cycle₁＝1,2,…,U₁) Taking ROI of the image at the intermediate time as a template frame Q_refPut it into the subset Q_exh,u1(ii) a In the same phase period of the next period, each frame of image is A × L in size, and the upper left corner is located at (a)₀+δa,l₀+ δ l) area is the detection frame Q_tar(ii) a Delta a and delta L are respectively changed within the range of-A/3 to A/3 and-L/3 to L/3, delta a and delta L are respectively axial displacement and transverse displacement, and Q is enabled to be_tarAnd Q_refThe two-dimensional cross correlation coefficient of (1) is maximum, and corresponding delta a and delta l are Q_tarThe displacement amounts Δ a and Δ l of; will Q_tarShifted in the axial and lateral directions by Δ a and Δ l, respectively, to obtain a displacement-corrected image Q'_tar(ii) a Comparing all maximum cross correlation coefficients, and comparing Q 'corresponding to the maximum value in the maximum cross correlation coefficients'_tarAs new Q_refPut into the subset Q_exh,u1To the end of (1); using new Q_refThe same operation is carried out on the next expiration period, and the operation is circulated until all expiration periods are calibrated and matched;

treating the inspiration period by the same method to obtain a subset Q_inh,u2(u₂＝1,2,…,U₂)。

Further, the specific implementation manner of the fusion step is as follows:

for image subset Q_uCalculating instantaneous displacement distribution between time-adjacent images by using a block matching algorithm on a thread-by-thread u basis; accumulating the instantaneous displacement distribution at different moments to obtain the spatial accumulated displacement D of each thread_u(ii) a By pair D_uCalculated along the axial direction to obtainAccumulated strain S_u；

Let each S_uSet of time coordinates t_u(U-1, 2, …, U), converting all t' s_uIs denoted by t_all(ii) a S obtained for each thread_uAnd (4) performing time-domain cubic spline interpolation to enable the strain result of each thread to cover a time coordinate set t_all；

Designing a denoising function based on the correlation among the threads of the accumulated strain, and denoising the strain distribution;

at t_allAt each time point in (2), the threads after noise reduction are accumulated with strain S'_uAveraging to obtain the final strain imaging result S_all。

Furthermore, the denoising function is designed based on the correlation between threads of the accumulated strain, and the specific implementation method for the strain distribution denoising is as follows:

at any time, calculate S_uCross correlation coefficients among threads; if the cross-correlation coefficient of a certain thread and other threads is lower than a threshold value gamma, the output strain of the thread at the current moment is eliminated, and the value range of gamma is 0.4-0.9;

at that moment, with the image S_uIs taken as the center, and the size is A₀×L₀Rectangular window R_u，A₀、L₀The axial and lateral dimensions of the rectangular window, in pixels, are the strain S of the other threads_v(v ≠ U ≠ 1,2, …) for the corresponding portion R_vThe average cross-correlation coefficient is calculated as follows,

wherein, COV (R)_u,R_v) Represents R_uAnd R_vCovariance of (a)_R1And σ_RvIs R_uAnd R_vStandard deviation of (2)

Denoising function q_u(i, j) is defined as:

wherein, alpha is attenuation coefficient, and the value range can be 10^-2～10²，p_thThe correlation coefficient threshold value is set, and the value range can be 0.4-0.9;

de-noising function q_u(i, j) and S_uMultiplying to obtain a denoised strain image S'_u。

The invention also provides a multithreading strain imaging device based on the living body ultrasonic image, which comprises the following components: the ultrasonic diagnostic system comprises a computer, a B ultrasonic instrument, a linear array probe and a strain generating device;

the B ultrasonic instrument is connected with the linear array probe and is used for acquiring ultrasonic images of the living body object;

strain generating means for creating a displacement and strain distribution in the tissue;

the computer is respectively connected with the B-ultrasonic instrument and the strain generating device, is used for controlling the B-ultrasonic instrument and the strain generating device, and is configured to realize the above multithreading strain imaging method based on the living body ultrasonic image.

3. Advantageous effects

Compared with the prior art, the invention has the beneficial effects that:

according to the multithreading strain imaging method based on the living body ultrasonic image, an ultrasonic image sequence is divided into a plurality of independent threads corresponding to different motion phases, image registration, displacement compensation and strain imaging are performed in a multithreading parallel mode, strain distribution under a plurality of motion phases with motion being suppressed can be obtained, and strain imaging depending on a single motion phase is avoided; denoising the strain based on the spatial correlation of the strain distribution of the multiple threads, and further improving the signal-to-noise ratio of the strain distribution and the robustness of the algorithm; the calculation of a plurality of threads in the strain imaging is independent and parallel to each other, and the strain imaging system has an efficient data processing and calculating structure; after strain denoising and fusion of multiple threads, the finally output strain has obviously improved precision and time density, and has extremely high application value for clinical real-time performance.

Drawings

FIG. 1 is a schematic flow chart of a B-mode ultrasound image-based in-vivo multithread strain imaging method;

FIG. 2 is a schematic view of the apparatus of example 2;

FIG. 3 shows p in example 2_inh(t) setting a threshold value of a second inspiration period, dividing motion phases and selecting images of all phases;

FIG. 4 is the distribution after registration of 4 threads in inhalation and exhalation states in example 2;

FIG. 5 is a cumulative thermal strain profile of a sequence of 4 thread images for the inspiratory phase and the expiratory phase of example 2;

FIG. 6 is a schematic diagram showing the spatial dependence of the thermal strain on each thread in example 2;

FIG. 7 is a diagram illustrating the cumulative thermal strain distribution of each thread after the denoising function is applied in example 2;

FIG. 8 is a strain-time curve of the maximum strain point after strain fusion for each thread in example 2;

reference numerals: 01. a computer; 02. a portable B ultrasonic instrument; 03. b ultrasonic imaging probe; 04. a microwave heater; 05. a microwave ablation needle.

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some embodiments of the present invention, but not all embodiments; moreover, the embodiments are not relatively independent, and can be combined with each other according to needs, so that a better effect is achieved. Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

For a further understanding of the invention, reference should be made to the following detailed description taken in conjunction with the accompanying drawings and examples.

Example 1

With reference to fig. 1, a multithreading strain imaging method based on a living body ultrasound image of the present invention includes the following steps:

an acquisition step, namely acquiring an ultrasonic image of a living object to obtain a digital ultrasonic image sequence I:

b-type ultrasonic imaging is carried out on biological tissues of a target area to be inspected aiming at a living object in a conventional physiological motion state, and data acquisition is carried out at equal time intervals to obtain a digital image sequence I;

specifically, the B ultrasonic instrument works in a real-time imaging mode, the position, the angle and the imaging parameters of the ultrasonic probe are adjusted, and the parameters of the probe and the system are fixed after the target area tissues are clearly visible in the image; by applying external force or heat, etc., displacement and strain distribution is formed in the tissue; b ultrasonic continuously and at equal time intervals collects digital images; the acquisition duration is not shorter than 2 breathing cycles.

The ultrasonic instrument for real-time imaging and B ultrasonic digital image acquisition can be a cart type B ultrasonic machine or a portable type B ultrasonic machine; the adopted ultrasonic imaging probe can be an abdominal probe, a craniocerebral probe, an intracavity probe and a heart probe according to the diagnosis part, can be a rectangular probe, an arc probe and a circular probe according to the geometrical shape, and can be a linear scanning probe, a phased array probe and a mechanical sector scanning probe according to the beam control mode; the parameters of the imaging system are selected from imaging depth, time gain compensation, focusing depth, emission frequency, linear density and the like; the beam transmitting logic and the corresponding beam forming algorithm of the B-ultrasonic instrument are selected according to the strain imaging algorithm, the motion rate of scattering points and the like and by integrating the requirements of the transverse resolution and the frame rate of the image, and the method comprises focused beam imaging, plane wave imaging, multi-beam imaging, synthetic aperture imaging and the like.

The sources for applying external force to the target area tissues comprise quasi-static pressure applied in vitro, elastic waves formed by low-frequency external driving mechanical force in the tissues, acoustic radiation force formed by the action of a focused ultrasonic sound field on the tissues, shear waves generated by acoustic pulses and transversely transmitted along beams, and the like; the heating mode comprises means such as radio frequency ablation, microwave ablation, laser ablation, focused ultrasound ablation and the like.

The stability of the living body state is kept in the process of collecting the image, and the breathing mode can be autonomous breathing, breathing supported by a breathing machine and the like; the heartbeat may be an autonomous heartbeat, a heartbeat supported by a cardiac pacemaker, or the like; the acquisition process at least acquires image sequences with 2 physiological motion cycles and more, so that effective spectrum analysis is facilitated; the data form of the B-ultrasonic digital image can be an ultrasonic Radio Frequency (RF) signal, an analytic signal of the RF signal after orthogonal demodulation or Hilbert transform, a gray scale signal obtained by enveloping and compressing the signal, and the like.

A characteristic extraction step: extracting the texture features of the digital ultrasonic image sequence I:

in this embodiment, correlation calculation between images is used as a feature analysis means for image texture information, and extracting texture features is not necessarily limited to extracting by correlation calculation, and may be performed based on information such as energy, contrast, entropy, and the like of images.

Respectively taking the image which is firstly acquired and covers a respiratory cycle in the digital ultrasonic image sequence I as a template frame I_ref,n(n is an image serial number), and the two-dimensional cross correlation coefficient calculation is carried out on the image serial number and the rest images to obtain a series of curves p of which the cross correlation coefficient changes along with time_n(t); specifically, N is 1,2, …, M frames in N frame images are taken as template frames I respectively_ref(N is the total frame number of the image sequence I, M is a positive integer and makes the M frame image cover at least one complete respiration cycle), and the rest N-1 frames are used as the detection frame I_tarCalculating I_refAnd I_tarObtaining M cross correlation coefficients-time curve p_n(t)。

The specific method of the two-dimensional correlation calculation is not limited to the normalized cross-correlation algorithm, and also includes a mean absolute difference algorithm, a sum of absolute differences algorithm, a sum of square errors algorithm, a mean sum of square errors algorithm, a sequential similarity test algorithm, non-normalized cross-correlation, Hybrid-sign correlation, polar-correlation or Meyr-Spies method and other known methods. Taking normalized cross-correlation as an example, let the templateThe matrix of frame pixels is I_refDetecting a frame pixel matrix as I_tarThen the cross correlation coefficient p is calculated as follows:

wherein, COV (I)_ref,I_tar) Is represented by_ref,I_tarThe covariance of (a) of (b),

is represented by_ref,I_tarStandard deviation of (2).

A cycle dividing step: determining a characteristic image of an expiration state and an inspiration state through the texture characteristics of the digital ultrasonic image sequence I, and dividing the digital ultrasonic image sequence I into an image sequence of a plurality of expiration periods and an image sequence of a plurality of inspiration periods:

within the interval of possible physiological movement frequencies, from p_n(t) frequency spectrum P_n(f) Characteristic image I for determining expiration and inspiration states_RE、I_RIAnd the respiratory frequency f_res(ii) a In particular, the frequency spectrum P_n(f) Is denoted as A_n，A_nThe template frame corresponding to the maximum value of (n-1, 2, …, M) is I_RECorresponding to p_n(t) is denoted by p_exh(t); to p_exh(t) after low-pass filtering, the detection frame corresponding to the minimum value in the first M points on the curve is I_RI，I_RIAs template frame correspondences_n(t) is denoted by p_inh(t)；p_exh(t) the frequency corresponding to the maximum of the spectral amplitude is f_res。

Time domain signal p_n(t) transforming to obtain a frequency spectrum P_n(f) To determine the breathing frequency requires high spectral resolution, methods of spectral refinement that may be employed include the refined fast fourier transform (zoomft), the chirped z-transform, and the Yip-Zoom transform, among others.

To I_RECorresponding correlation coefficient-time curve p_exh(t) low-pass filtering, the time corresponding to all peak values of the filtered curveTaking the moment as an inspiration period boundary and the corresponding moment of the valley value as an expiration period boundary, and dividing all the images into different periods of expiration and inspiration states according to the moment; in particular, the FIR low-pass filter pair p is designed_exh(t) Filtering to obtain a smoothed curve p'_exh(t) the filter cut-off frequency is selected at f_res～3f_resTo (c) to (d); p'_exh(t) an inspiration period is formed between every two adjacent peaks, and an expiration period is formed between the valleys; additionally, there is an incomplete period from the starting time to the first extreme, and if the extreme is a peak (valley), there is an incomplete inspiration (expiration) period; the last extreme to the end is an incomplete cycle, and if the extreme is a peak (valley), it is an incomplete inspiratory (expiratory) cycle.

From p_exh(t) smoothing to give p'_exhThe method of (t) is not limited to FIR low pass filters, but other smoothing filters may be used, including mean filtering, median filtering, local regression, Savitzky-Golay filtering, and the like.

For low pass filter post curve p'_exh(t) poor smoothness may result in erroneous peak and valley selection, which may affect accurate cycle division. Can be obtained by imaging the exhalation characteristics I_RESet as an anchor point at a breathing frequency f_resDesigning the search step Δ t, e.g. 1/f_res>Δt>0.5/f_resIn the direction of increasing and decreasing time from anchor point, at p'_exh(t) alternately searching for valleys and peaks.

Setting a threshold value: setting a expiration textural feature threshold, and if the textural feature of the image in the expiration period is greater than the expiration textural feature threshold, placing the image in an expiration image sequence I_exhDividing the image sequence of each expiratory cycle into U phases₁A phase period, U₁For the number of divided phase periods, a multithreading subset I is obtained_exh,u1(u₁＝1,2,…,U₁) (ii) a Setting an inspiration textural feature threshold, and if the textural feature of the image in the inspiration period is larger than the inspiration textural feature threshold, placing the image in an inspiration image sequence I_inhFor each inspirationThe periodic image sequence is divided into U phases by dividing the motion phases₂A phase period, U₂For the number of divided phase periods, a multithreading subset I is obtained_inh,u2(u₂＝1,2,…,U₂)：

At p_exh(t) in each expiratory cycle, setting the upper quartile of all cross-correlation coefficients as the expiratory texture feature threshold, for p_exh(t) points above the expiratory texture feature threshold, and placing the corresponding images into the expiratory image sequence I_exh；

Dividing the image sequence of each expiration period into U phases₁A phase period, U₁For the number of divided phase periods, a multithreading subset I is obtained_exh,u1(u₁＝1,2,…,U₁)；

At p_inh(t) in each inspiration period, setting the upper quartile of all cross-correlation coefficients as the inspiration texture feature threshold, for p_inh(t) points above the threshold of the inspiratory texture feature, and placing the corresponding images in the inspiratory image sequence I_inh；

Dividing the image sequence of each inspiration period into U₂A phase period, U₂For the number of divided phase periods, a multithreading subset I is obtained_inh,u2(u₂＝1,2,…,U₂)。

The method for dividing the motion phase of the image sequence of each expiration period is to divide or set the division of the amplitude interval of the correlation coefficient at equal time intervals. Generally, the time of inspiration and expiration states of the living body are not equal, so the number of phases U of the two states₁、U₂Nor are they necessarily equal.

The texture feature threshold is selected from an upper quartile of cross correlation coefficients from small to large in the period, and the selection of the texture feature threshold is related to the overall correlation level of the image sequence and is not necessarily limited to the upper quartile.

The threshold setting is for p_exh(t) and p_inh(t) incomplete length of cycles when the incomplete length of exhalationOr the length of the inspiration period is less than 0.9/f_resThen the threshold for that cycle is taken to be the threshold for the adjacent full cycle.

A registration step: determining the region of interest by taking the position of the target region as the center, for I_exh,u1、I_inh,u2To the region of interest, displacement calibration and image registration are performed to obtain registered subsets Q_exh,u1、Q_inh,u2；Q_exh,u1And Q_inh,u2Is denoted as Q_u(U-1, 2, …, U) and U is Q_uNumber of threads of (U) ═ U₁+U₂：

The specific process comprises the following steps: the position of the target area is taken as the center, the size is determined to be A (axial direction) multiplied by L (transverse direction), and the coordinate of the upper left corner is determined to be (a)₀,l₀) ROI of interest.

During each phase (u) of the first expiratory cycle₁＝1,2,…,U₁) Taking ROI of the image at the intermediate time as a template frame Q_refPut it into the subset Q_exh,u1(ii) a In the same phase period of the next period, each frame of image is A × L in size, and the upper left corner is located at (a)₀+δa,l₀+ δ l) area is the detection frame Q_tar(ii) a Delta a and delta L are respectively changed within the range of-A/3 to A/3 and-L/3 to L/3, delta a and delta L are respectively axial displacement and transverse displacement, and Q is enabled to be_tarAnd Q_refThe two-dimensional cross correlation coefficient of (1) is maximum, and corresponding delta a and delta l are Q_tarThe displacement amounts Δ a and Δ l of; will Q_tarShifted in the axial and lateral directions by Δ a and Δ l, respectively, to obtain a displacement-corrected image Q'_tar(ii) a Comparing all maximum cross correlation coefficients, and comparing Q 'corresponding to the maximum value in the maximum cross correlation coefficients'_tarAs new Q_refPut into the subset Q_exh,u1To the end of (c).

Using new Q_refThe same is done for the next expiratory cycle, and so on until all expiratory cycles are calibrated and matched. Treating the inspiration period by the same method to obtain a subset Q_inh,u2(u₂＝1,2,…,U₂). All Q_exh,u1And Q_inh,u2The collection of (total U threads) is marked as Q_u。

A fusion step:parallel computation Q_uPerforming time domain interpolation on the accumulated strain image sequence of each thread, denoising based on strain correlation, fusing the accumulated strain of each thread after processing to obtain a final strain imaging result:

for image subset Q_uCalculating instantaneous displacement distribution between time-adjacent images by using a block matching algorithm on a thread-by-thread u basis; accumulating the instantaneous displacement distribution at different moments to obtain the spatial accumulated displacement D of each thread_u(ii) a By pair D_uThe cumulative strain S is calculated by the difference along the axial direction (the direction of the ultrasonic beam)_u。

Performing time domain interpolation on the accumulated strain image sequence of each thread to enable all threads to have the same time coordinate; specifically, let each S_uSet of time coordinates t_u(U-1, 2, …, U), converting all t' s_uIs denoted by t_all(ii) a S obtained for each thread_uAnd (4) performing time-domain cubic spline interpolation to enable the strain result of each thread to cover a time coordinate set t_all。

Designing a denoising function based on the correlation between threads of accumulated strain, wherein the specific process of denoising the strain distribution comprises the following steps:

at any time, calculate S_uCross correlation coefficients among threads; if the cross-correlation coefficient of a certain thread and other threads is lower than a threshold value gamma, the output strain of the thread at the current moment is eliminated, and the value range of gamma is 0.4-0.9;

at that moment, with the image S_uIs taken as the center, and the size is A₀×L₀Rectangular window R_u，A₀、L₀The axial and lateral dimensions of the rectangular window, in pixels, are the strain S of the other threads_v(v ≠ U ≠ 1,2, …) for the corresponding portion R_vThe average cross-correlation coefficient is calculated as follows,

wherein, COV (R)_u,R_v) Represents R_uAnd R_vCovariance of (a)_R1And σ_RvIs R_uAnd R_vStandard deviation of (2). Denoising function q_u(i, j) is defined as:

wherein, alpha is attenuation coefficient, and the value range can be 10^-2～10²，p_thThe correlation coefficient threshold value is a value range of 0.4-0.9.

De-noising function q_u(i, j) and S_uMultiplying to obtain a denoised strain image S'_u。

Accumulating strain S 'for each thread subjected to noise reduction'_uFusing to obtain final strain imaging result S_allIn the case of a liquid crystal display device, in particular,

fusing the strain results of all threads to obtain a final strain imaging result S_allThe specific process comprises the following steps: at t_allAt each time point in (2), the cumulative strain distribution obtained by each thread is averaged.

It is worth explaining that, according to the multi-thread strain imaging method based on the living body ultrasonic image, the more phases/threads are selected, the higher the accuracy, the time density and the stronger the real-time performance of the finally output strain result are, and the robustness of the strain imaging algorithm is stronger by combining the space correlation denoising of the thread results. Theoretically, the algorithm can cover all phases of respiratory motion, but in practical application, the algorithm is limited by a limited ultrasonic image frame rate, and physiological motion is not a complete ideal period, so that motion phases are divided in an interval where respiratory motion is relatively stable, and the number of selectable motion phases is not infinite; in addition, the method of the invention is essentially to process the live images which are arranged in sequence, so that the method is definitely applicable to two-dimensional live image sequences with similar characteristics obtained by other physical means besides B-ultrasonic.

It should be further noted that the invention is a kind of image data driven door-like motion suppression method, under the condition of not depending on the external equipment, have realized the state recognition of physiological motion, physiological motion frequency discrimination, physiological motion cycle division, have greatly simplified the monitoring equipment; the image sequences of all threads of the invention are independent in the links with relatively complex calculation, such as image registration, displacement correction and strain calculation, so that the parallel programming structure is adopted, and the calculation efficiency of the strain imaging algorithm can be obviously improved. After each thread completes independent strain calculation, the strain output which is dense in time is obtained through fusion, and the method has high application value for real-time monitoring in clinical treatment.

Example 2

In this embodiment, the method in embodiment 1 is adopted, the microwave ablation needle is used to heat the visceral fat of the living pig and acquire the B-ultrasonic digital image signal, and the thermal strain of the living pig is calculated by the B-ultrasonic image-based multithread strain imaging method, which includes the following steps:

a multithread strain imaging device based on living body ultrasonic images is shown in figure 2, a computer 01 controls digital ultrasonic image acquisition of a portable B ultrasonic instrument 02, the digital sampling rate of the B ultrasonic instrument is 40MHz, an adopted imaging probe is a 128-array-element linear array probe 03, and the center frequency of the imaging probe is 10.5 MHz. Under the guidance of a real-time B ultrasonic image presented by a B ultrasonic instrument, adjusting the imaging position and angle of an ultrasonic probe and the parameters of an imaging system to clearly present target area tissues in the center of the B ultrasonic image; the strain generating device is used to form displacement and strain distribution in the tissue, and the strain generating device forms strain by external force action, heating, or the like, and the specific mode is not limited. In the embodiment, displacement and strain distribution are formed by adopting a heating mode, a microwave heater 04 is connected with a microwave ablation needle 05 to serve as a heating source, and the microwave ablation needle is inserted into adipose tissues of a target region in a direction normal to an ultrasonic imaging plane along a wound of the epidermis of a live pig under the guidance of a B ultrasonic image, so that the heating region is in the imaging plane; the imaging depth of the B-mode ultrasound machine was set to 4cm, and a digital ultrasound image sequence of 2048 (axial) × 128 (transverse) pixels per frame was acquired at a frame rate of 50Hz for a total of 20 seconds from the start of heating.

Heating process 20 seconds, collecting N-1000 frames of images, making M-250 to ensure covering a complete respiration period, adopting normalized cross-correlation to calculate 1-M frames of B-ultrasonic images as template frame I_refAnd the rest N-1 frame images are taken as detection frames I_tarCross correlation coefficient p of_n(t)。

In the target frequency band f₁＝0.1Hz，f₂In the interval 5Hz, p is obtained_n(t) refined frequency spectrum P_n(f) According to P_n(f) The amplitude relation can determine the expiration characteristic image I_REFor frame 92, inspiratory phase feature image I_RIFor frame 396, respiratory rate f_res＝0.3325Hz。

To I_RECorresponding correlation coefficient-time curve p_exh(t) low-pass filtering with FIR low-pass filter with normalized cut-off frequency of 0.02 and order of 81, and compensating group delay to obtain p'_exh(t) of (d). From p'_exh(t) the peak value and the valley value of the (t) are obtained to obtain the periodic boundary t of the expiratory state _exh0,2.46,5.28,7.94,10.98,13.7,16.58,19.46,20s, and cycle boundary t of inspiratory state_inhThe division results in 8 expiration periods and 8 inspiration periods, namely 0,1.04,3.86,6.6,9.46,12.38,15.18,18.14 and 20 s.

At p_exh(t) in each expiration cycle, setting the upper quartile of all cross-correlation coefficients as a threshold, and placing the corresponding template frame larger than the threshold into I_exh(ii) a Similarly, at p_inh(t) in each inspiration period, setting the upper quartile of all cross-correlation coefficients as a threshold, and placing the corresponding template frame larger than the threshold into I_inh(ii) a Each expiratory cycle is divided equally into U₁＝U ₂4 different phases for 4 phase periods, resulting in a multi-thread subset I_exh,u1(u₁＝1,2,…,U₁)、I_inh,u2(u₂＝1,2,…,U₂). FIG. 3 shows a graph at p_inh(t) threshold setting for the second inspiratory cycle, motion phase division and image selection for each phase.

Taking the target area position of the B-ultrasonic image as the center,a region of interest of size 561 pixels x 35 pixels is determined. To I_exh,u1、I_inh,u2To the region of interest, displacement calibration and image registration are performed to obtain registered subsets Q_exh,u1、Q_inh,u2. All U being U₁+U₂The registered subset corresponding to each thread is named Q_u(u＝1,2,…,U₁+U₂). As shown in FIG. 4, FIGS. 4(a) -4(d) show the 4-thread registered subset Q in the inspiratory state_inh,u2At p_inh(t) distribution, FIGS. 4(e) -4(h) show the 4-thread registered subset Q in the expiratory state_exh,u1At p_exh(t) distribution over (t).

Image subset Q_uCalculating the displacements D in parallel_uStrain S_uThe cumulative thermal strain results when t is 20s are shown in fig. 5, where fig. 5(a) to 5(d) show the cumulative thermal strain profiles for the respective threads of the inhalation cycle, and fig. 5(e) to 5(h) show the cumulative thermal strain profiles for the respective threads of the exhalation cycle. Each thread develops a cumulative thermal strain of about 0.2 and has a similar thermal strain profile, with strain noise in the corners of the regions shown by the individual threads.

After the thermal strain time domain interpolation of each thread, the spatial correlation among strains is calculated to obtain the spatial correlation distribution

As shown in fig. 6, the region with high gray value represents high correlation in the space, and the region with low gray value represents low correlation. De-noising function q is designed based on the distribution_u(i,j)，p_th0.95, 5 and the strain distribution S of each thread_uMultiplication. As shown in fig. 7, fig. 7(a) -7(d) are cumulative thermal strain distributions after denoising for each thread in an inspiration period, and fig. 7(e) -7(h) are cumulative thermal strain distributions after denoising for each thread in an expiration period, and it can be seen that the strain noise outside the heat source is effectively suppressed.

Finally, accumulating strain S 'for each thread after noise reduction'_uFusing to obtain final strain imaging result S_allFIG. 8 shows the post-fusion thermal strain S_allThermal strain-time curve at the position of the maximum, scatter point isThe strain point of the output, the strain average output time interval is 0.35 second, the maximum output time interval is 0.76 second, and the time density is very high.

The invention has been described in detail hereinabove with reference to specific exemplary embodiments thereof. It will, however, be understood that various modifications and changes may be made without departing from the scope of the invention as defined in the appended claims. The detailed description and drawings are to be regarded as illustrative rather than restrictive, and any such modifications and variations are intended to be included within the scope of the present invention as described herein. Furthermore, the background is intended to be illustrative of the state of the art as developed and the meaning of the present technology and is not intended to limit the scope of the invention or the application and field of application of the invention.

Claims (10)

1. A multithreading strain imaging method based on a living body ultrasonic image is characterized by comprising the following steps:

collecting an ultrasonic image of a living object to obtain a digital ultrasonic image sequence I;

a characteristic extraction step: extracting the texture features of the digital ultrasonic image sequence I;

a cycle dividing step: determining characteristic images of the expiration and inspiration states through the texture characteristics of the digital ultrasonic image sequence I, and dividing the digital ultrasonic image sequence I into a plurality of image sequences of expiration periods and a plurality of image sequences of inspiration periods;

setting a threshold value: setting an expiration texture feature threshold;

if the texture feature of the image on the expiration period is larger than the expiration texture feature threshold value, the image is placed into an expiration image sequence I_exhDividing the image sequence of each expiratory cycle into U phases₁A phase period, U₁For the number of divided phase periods, a multithreading subset I is obtained_exh,u1(u₁＝1,2,…,U₁) (ii) a Setting a suction texture feature threshold;

if the texture of the image over an inspiratory period is greater than an inspiratory texture threshold, the image is placed inInspiratory image sequence I_inhDividing the image sequence of each inspiration cycle into U phases₂A phase period, U₂For the number of divided phase periods, a multithreading subset I is obtained_inh,u2(u₂＝1,2,…,U₂)；

A registration step: determining the region of interest by taking the position of the target region as the center, for I_exh,u1、I_inh,u2To the region of interest, displacement calibration and image registration are performed to obtain registered subsets Q_exh,u1、Q_inh,u2；Q_exh,u1And Q_inh,u2Is denoted as Q_u(U-1, 2, …, U) and U is Q_uNumber of threads of (U) ═ U₁+U₂；

A fusion step: parallel computation Q_uThe time domain interpolation is carried out on the accumulated strain image sequence of each thread, denoising is carried out based on strain correlation, and the accumulated strain of each thread after processing is fused to obtain a final strain imaging result.

2. The method of claim 1, wherein the texture features in the feature extraction step are correlations between images, and the specific extraction method comprises:

taking N-1, 2, …, M frames of N frames of image of digital ultrasonic image sequence I as template frame I_refN is the total frame number of the image sequence I, M is a positive integer, and the M frame images at least cover a complete respiration period, and the rest N-1 frames are used as detection frames I_tarCalculating I_refAnd I_tarObtaining M cross correlation coefficients-time curve p_n(t)。

3. The method for multi-thread strain imaging based on living body ultrasound images as claimed in claim 1 or 2, wherein the step of dividing the cycle is implemented by:

within the interval of possible physiological movement frequencies, from p_n(t) ofFrequency spectrum P_n(f) Characteristic image I for determining expiration and inspiration states_RE、I_RIAnd the respiratory frequency f_resThe calculation method comprises the following steps: will spectrum P_n(f) Is denoted as A_n，A_nThe template frame corresponding to the maximum value of (n-1, 2, …, M) is I_RECorresponding to p_n(t) is denoted by p_exh(t) for p_exh(t) after low-pass filtering, the detection frame corresponding to the minimum value in the first M points on the curve is I_RI，I_RIAs template frame correspondences_n(t) is denoted by p_inh(t)，p_exh(t) the frequency corresponding to the maximum of the spectral amplitude is f_res(ii) a To I_RECorresponding correlation coefficient-time curve p_exh(t) filtering to obtain a smooth curve p'_exhAnd (t) taking the corresponding time of all the peaks as the boundary of the inspiration period and the corresponding time of the valleys as the boundary of the expiration period, thereby dividing the digital ultrasonic image sequence I into a plurality of image sequences of expiration periods and a plurality of image sequences of inspiration periods.

4. The live ultrasound image-based multi-threaded strain imaging method of claim 3, wherein p 'is selected'_exhThe method of the peak and the valley of (t) is: the exhalation characteristic image I_RESet as an anchor point at a breathing frequency f_resDesign search step Δ t, 1/f_res>Δt>0.5/f_resIn the direction of increasing and decreasing time from anchor point, at p'_exh(t) alternately searching for valleys and peaks.

5. The method of claim 1, wherein the dividing of the motion phase of the image sequence for each expiratory cycle in the thresholding step is performed by dividing or setting the division of the correlation coefficient amplitude interval at equal time intervals.

6. The method of claim 1, wherein the threshold setting step is implemented by:

at p_exh(t) in each expiratory cycle, setting the upper quartile of all cross-correlation coefficients as the expiratory texture feature threshold, for p_exh(t) points above the expiratory texture feature threshold, and placing the corresponding images into the expiratory image sequence I_exh；

Dividing the image sequence of each expiration period into U phases₁A phase period, U₁For the number of divided phase periods, a multithreading subset I is obtained_exh,u1(u₁＝1,2,…,U₁)；

At p_inh(t) in each inspiration period, setting the upper quartile of all cross-correlation coefficients as the inspiration texture feature threshold, for p_inh(t) points above the threshold of the inspiratory texture feature, and placing the corresponding images in the inspiratory image sequence I_inh；

Dividing the image sequence of each inspiration period into U₂A phase period, U₂For the number of divided phase periods, a multithreading subset I is obtained_inh,u2(u₂＝1,2,…,U₂)。

7. The method of claim 1, wherein the registering step is implemented as follows:

the position of the target area is taken as the center, the size is determined to be AxL, and the coordinate of the upper left corner is determined to be (a)₀,l₀) Wherein a is the axial length and L is the transverse length;

during each phase (u) of the first expiratory cycle₁＝1,2,…,U₁) Taking ROI of the image at the intermediate time as a template frame Q_refPut it into the subset Q_exh,u1(ii) a In the same phase period of the next period, each frame of image is A × L in size, and the upper left corner is located at (a)₀+δa,l₀+ δ l) area is the detection frame Q_tar(ii) a Delta a and delta L respectively varying in the range-A/3 to A/3 and-L/3 to L/3, delta a and delta L respectively being axial and transverseTo the displacement amount of Q_tarAnd Q_refThe two-dimensional cross correlation coefficient of (1) is maximum, and corresponding delta a and delta l are Q_tarThe displacement amounts Δ a and Δ l of; will Q_tarShifted in the axial and lateral directions by Δ a and Δ l, respectively, to obtain a displacement-corrected image Q'_tar(ii) a Comparing all maximum cross correlation coefficients, and comparing Q 'corresponding to the maximum value in the maximum cross correlation coefficients'_tarAs new Q_refPut into the subset Q_exh,u1To the end of (1); using new Q_refThe same operation is carried out on the next expiration period, and the operation is circulated until all expiration periods are calibrated and matched;

during each phase (u) of the first inspiratory cycle₂＝1,2,…,U₂) Taking ROI of the image at the intermediate time as a template frame Q_refPut it into the subset Q_inh,u2(ii) a In the same phase period of the next period, each frame of image is A × L in size, and the upper left corner is located at (a)₀+δa,l₀+ δ l) area is the detection frame Q_tar(ii) a Delta a and delta L are respectively changed within the range of-A/3 to A/3 and-L/3 to L/3, delta a and delta L are respectively axial displacement and transverse displacement, and Q is enabled to be_tarAnd Q_refThe two-dimensional cross correlation coefficient of (1) is maximum, and corresponding delta a and delta l are Q_tarThe displacement amounts Δ a and Δ l of; will Q_tarShifted in the axial and lateral directions by Δ a and Δ l, respectively, to obtain a displacement-corrected image Q'_tar(ii) a Comparing all maximum cross correlation coefficients, and comparing Q 'corresponding to the maximum value in the maximum cross correlation coefficients'_tarAs new Q_refPut into the subset Q_inh,u2To the end of (1); using new Q_refThe same is done for the next inspiratory cycle, and so on until all inspiratory cycles are calibrated and matched.

8. The method of claim 1 or 7, wherein the fusion step is implemented as follows:

for image subset Q_uCalculating instantaneous displacement distribution between time-adjacent images by using a block matching algorithm on a thread-by-thread u basis; accumulating the instantaneous displacement distribution at different moments to obtain the spatial accumulated displacement D of each thread_u(ii) a By pair D_uThe accumulated strain S is obtained by calculating the difference along the axial direction_u；

Let each S_uSet of time coordinates t_u(U-1, 2, …, U), converting all t' s_uIs denoted by t_all(ii) a S obtained for each thread_uAnd (4) performing time-domain cubic spline interpolation to enable the strain result of each thread to cover a time coordinate set t_all；

Designing a denoising function based on the correlation among the threads of the accumulated strain, and denoising the strain distribution;

at t_allAt each time point in (2), the threads after noise reduction are accumulated with strain S'_uAveraging to obtain the final strain imaging result S_all。

9. The method as claimed in claim 8, wherein the method for multi-threaded strain imaging based on living body ultrasound images comprises the following steps:

at any time, calculate S_uCross correlation coefficients among threads; if the cross-correlation coefficient of a certain thread and other threads is lower than a threshold value gamma, the output strain of the thread at the current moment is eliminated, and the value range of gamma is 0.4-0.9;

at said arbitrary time, with the image S_uIs taken as the center, and the size is A₀×L₀Rectangular window R_u，A₀、L₀The axial and lateral dimensions of the rectangular window, in pixels, are the strain S of the other threads_v(v ≠ U ≠ 1,2, …) for the corresponding portion R_vThe average cross-correlation coefficient is calculated as follows,

wherein, COV (R)_u,R_v) Represents R_uAnd R_vCovariance of (a)_R1And σ_RvIs R_uAnd R_vStandard deviation of (2)

Denoising function q_u(i, j) is defined as:

wherein alpha is attenuation coefficient and has a value range of 10^-2～10²，p_thThe correlation coefficient threshold value is a value range of 0.4-0.9;

de-noising function q_u(i, j) and S_uMultiplying to obtain a denoised strain image S'_u。

10. A multi-threaded strain imaging apparatus based on an ultrasound image of a living body, comprising: the ultrasonic diagnostic system comprises a computer, a B ultrasonic instrument, a linear array probe and a strain generating device;

the B ultrasonic instrument is connected with the linear array probe and is used for acquiring ultrasonic images of the living body object;

strain generating means for creating a displacement and strain distribution in the tissue;

a computer connected to the B-ultrasonic apparatus and the strain generating device, respectively, for controlling the B-ultrasonic apparatus and the strain generating device, configured to implement a method of multi-threaded strain imaging based on live ultrasound images as claimed in any one of claims 1 to 9.

Images