Methods for performing high spatiotemporal resolution ultrasound imaging of microvessels

Abstract

To perform high spatial and temporal resolution ultrasound imaging of microvessels in a subject.SOLUTION: Ultrasound data is acquired from a region-of-interest in a subject who has been administered a microbubble contrast agent. The ultrasound data is acquired while the microbubbles are moving through, or otherwise present in, the region-of-interest. The region-of-interest may include, for instance, microvessels or other microvasculature in the subject. By imaging microbubbles, a cross-correlation map between each microbubble image and a point spread function of the system can be generated. Accumulation of power-based cross-correlation maps may then be used to generate a high-resolution high-contrast image of the microvasculature.SELECTED DRAWING: Figure 1

Description

Technical Field

[0001] <Cross - reference to Related Applications> This application claims the benefit of U.S. Provisional Patent Application No. 63 / 039,549, filed on June 16, 2020, with the invention name "Methods for High Spatial and Temporal Resolution Ultrasound Imaging of Microvessels", and the entire content of the same is incorporated herein by reference.

[0002] <Statement Regarding Federally Sponsored Research> This invention was made with government support under Grant No. NS111039 awarded by the National Institutes of Health. The government has certain rights in this invention.

Background Art

[0003] Ultrasound contrast microbubbles are commonly used in clinical settings for ultrasound imaging enhancement. The typical size of microbubbles ranges from 1 to 5 μm, which is significantly smaller than the wavelength of ultrasound, typically on the order of 100 - 800 μm. As a result, individual microbubbles appear as point sources of a size comparable to the ultrasound wavelength, creating blurring. Also, microbubbles smaller than half of the ultrasound wavelength cannot be resolved individually, resulting in what is known as the diffraction limit of ultrasound.

[0004] Contrast-enhanced ultrasound (CEUS) imaging is now widely used in clinical practice to provide highly sensitive blood flow and perfusion imaging in both superficial and deep tissues with very high temporal resolution (real time). While CEUS can provide quantitative assessment of relative blood volume and tissue perfusion, it still has limitations in spatial resolution, and at that spatial resolution, it is almost impossible to distinguish microvascular structures. For example, many pathological changes in tissue, such as tumors, occur at the microvascular level, so imaging of vascular structures and hemodynamic properties of blood flow at the microvascular or capillary level (e.g., less than 100 μm) is desired in clinical practice.

[0005] Ultra-high-resolution ultrasound technology is already known, enabling imaging of microvascular structures beyond the ultrasound diffraction limit. This imaging is based on the pinpoint localization of individual microbubbles (MBs) and tracking of their movement. While spatial resolution is improved, a major problem remains: insufficient temporal resolution when microvascular images must be generated with data acquisition in tens of seconds. This significantly hinders the practical application of this technology in clinical settings. Clinically relevant acquisition times can be on the order of one or two seconds. In such timeframes, ultra-high-resolution technology may not be able to provide a completely comprehensive vascular lumen image, and attempting to force such a result can lead to poor visualization of microvessels. Furthermore, ultra-high-resolution methods typically involve complex localization and MB pairing. This relies on tracking algorithms, which incur significant computational costs.

[0006] Therefore, there is still a demand for high-resolution ultrasound microvascular imaging that can simultaneously achieve both high spatial and temporal resolution. There is also a demand for image generation in a very short time, comparable to clinical CEUS procedures. [Overview of the project]

[0007] This disclosure relates to the use of an ultrasound system to perform high spatiotemporal resolution imaging of microvessels. This method addresses the aforementioned shortcomings. Ultrasound data can be accessed from an image storage system, or ultrasound data can be acquired in other ways from the region of interest of a subject administered a microbubble contrast agent. Ultrasound data is acquired while the microbubbles are moving through the region of interest, or in other ways while they are present within the region of interest. The region of interest may include, for example, the subject's microvessels or other microvascular systems. High-resolution, high-contrast images of microvascular systems can be generated using a cross-correlation map between the microbubble images and the system's point spread function. Velocity estimation can also be performed from the microbubble trajectories obtained across multiple frames of the cross-correlation map.

[0008] One embodiment provides a method for performing high spatiotemporal resolution imaging of microvessels using an ultrasound system. The method includes the step of accessing, using a computer system, ultrasound data acquired using the ultrasound system from a region of interest of a subject in which a microbubble contrast agent was present when the ultrasound data was acquired. The method also includes the step of generating microbubble signal data by separating the microbubble signals in the ultrasound data from other signals in the ultrasound data using the computer system. The method further includes the steps of determining the point spread function (PSF) of the ultrasound system used to acquire the ultrasound data, and generating a cross-correlation map between the microbubble signal data and the point spread function of the ultrasound system. High spatiotemporal resolution microvessel images can be generated based at least partially on the generated cross-correlation map.

[0009] In one embodiment, an ultrasound system is provided for high spatiotemporal resolution imaging of microvessels. The system includes a computer system, which is configured to access the ultrasound data acquired using the ultrasound system from the region of interest of a subject in which a microbubble contrast agent was present when the ultrasound data was acquired. The computer system is further configured to generate microbubble signal data by separating the microbubble signals in the ultrasound data from other signals in the ultrasound data. The computer system is also configured to determine the point spread function (PSF) of the ultrasound system used to acquire the ultrasound data and to generate a cross-correlation map between the microbubble signal data and the point spread function of the ultrasound system. Based at least partially on the generated cross-correlation map, the computer system can generate high spatiotemporal resolution microvessel images.

[0010] One embodiment provides a method for performing high spatiotemporal resolution imaging of microvessels using an ultrasound system. The method includes the step of accessing, using a computer system, ultrasound data acquired using the ultrasound system from a region of interest of a subject in which a microbubble contrast agent was present when the ultrasound data was acquired. The method also includes the step of generating microbubble signal data by separating the microbubble signals in the ultrasound data from other signals in the ultrasound data using the computer system, and the step of separating the individual microbubble trajectories in the microbubble signal data. The method further includes the step of determining the trajectory parameters of the separated microbubble trajectories, and the step of estimating the microbubble velocity based on the determined trajectory parameters of the separated microbubble trajectories. High spatiotemporal resolution microvessel images can be generated, at least partially based on the estimated microbubble velocity.

[0011] The above and other aspects and advantages of the present disclosure will become apparent from the following description. The following description refers to the accompanying drawings, which constitute part of the description and show preferred embodiments for illustrative purposes. However, these embodiments do not necessarily represent the entire scope of the present invention. Therefore, when interpreting the scope of this invention, the claims should be taken into consideration. [Brief explanation of the drawing]

[0012] [Figure 1] This is a flowchart outlining the steps of an example method for generating ultra-high-resolution images of microvascular systems. [Figure 2] This figure shows an example of a microvessel that moves across multiple consecutive time frames due to tissue movement. [Figure 3] This flowchart outlines the steps of an example method for generating high-resolution microvascular images based on a power-based cross-correlation method. [Figure 4A] This is a non-restrictive example of a cross-correlation calculation between a sample data set and a point spread function for a one-dimensional Gaussian distribution with an Nth power. [Figure 4B] Figure 4A shows a non-restrictive example of the relationship between the power of the cross-correlation and FWHM. [Figure 5] This is an unrestricted example of the cross-correlation between one example of data and another N-power point spread function, calculated using multiple different N-power 1-dimensional Gaussian distributions. [Figure 6A] This is an image of the original MB ultrasound image data for a non-restrictive example. [Figure 6B] Figure 6A shows a non-restrictive example of a first-square-based cross-correlation map. [Figure 6C] Figure 6A shows a non-restrictive example of a cross-correlation map based on the fourth power. [Figure 6D] Figure 6A shows an image of the cross-correlation map based on the 8th power, as a non-restrictive example. [Figure 7A] This is a non-limiting example of an image when cross-correlation threshold sharpening is not performed. [Figure 7B] This is a non-limiting example of an image when the threshold sharpening process of cross-correlation is performed with the value set to 0.25. [Figure 7C] This is a non-limiting example of an image when the threshold sharpening process of cross-correlation is performed with the value set to 0.5. [Figure 7D] This is a non-limiting example of an image when the threshold sharpening process of cross-correlation is performed with the value set to 0.75. [Figure 8] This is a diagram showing the original MB data and cross-correlation map of a non-limiting example over multiple frames, and the high-definition blood flow image obtained by integration therefrom. [Figure 9] This is a diagram showing the original image and microvascular image of a non-limiting example processed with pancreatic tumor and liver tissues. [Figure 10] This is a diagram showing the original image and the thresholded image of a non-limiting example. [Figure 11] (A) is a diagram showing a non-limiting example of an MB trajectory in a three-dimensional spacetime having two spatial dimensions and one temporal dimension. (B) is a detailed view of a part of (A). (C) is a diagram showing a non-limiting example of the separated MB trajectory in the three-dimensional spacetime. (D) is a diagram showing a non-limiting example of the separated MB trajectory projected onto the axial-time plane. (E) is a diagram showing a non-limiting example of the separated MB trajectory projected onto the lateral-time plane. [Figure 11F] This is a flowchart describing each step of an example method for generating a high-definition hemodynamic microvascular image based on the MB velocity estimation method. [Figure 12A] This is the original MB image of a non-limiting example of a spatial dependence threshold discrimination method based on noise floor estimation. [Figure 12B] This is a non-limiting example of a measured noise floor image of a non-limiting example of a spatial dependence threshold discrimination method based on noise floor estimation. [Figure 12C] This is a non-limiting example of a noise floor image of a non-limiting example of a spatial dependence threshold discrimination method based on noise floor estimation. <0OO0101>This is a noise floor image of another unrestricted example of a spatially dependent thresholding method based on noise floor estimation. [Figure 12E] This is a non-specific example image before residual noise suppression processing. [Figure 12F] This is a non-specific example image after residual noise suppression processing. [Figure 13] This is a block diagram of an example ultrasonic system that can implement the method described in the present disclosure. [Figure 14] This is a block diagram of an example computer system capable of implementing the embodiments described in the present disclosure. [Modes for carrying out the invention]

[0013] The following describes a system and method for performing high spatiotemporal resolution ultrasound imaging of a subject's microvessels. High-resolution microvessel imaging in a subject can be used to image microvessels with high resolution and a high contrast-to-noise ratio (CNR). Conventional ultrasound imaging and CEUS imaging generally display images as signal intensity. In this disclosure, a high spatiotemporal resolution ultrasound imaging technique can generate images with significantly improved spatial resolution and CNR using correlation maps. Ultrasound data is acquired from a region of interest in a subject administered a microbubble contrast agent. The ultrasound data is acquired while the microbubbles are moving through the region of interest or while the microbubbles are present within the region of interest in other ways. The region of interest may include, for example, the subject's microvessels or other microvascular systems. By imaging the microbubbles (MBs), a cross-correlation map can be generated between each MB image and the system's point spread function (PSF). Subsequently, by integrating the power-based cross-correlation maps over a certain number of frames, high-resolution, high-contrast images of microvessels can be generated, as described in this disclosure.

[0014] Referring here to Figure 1, a flowchart is shown illustrating each step of an example method using an ultrasound system to generate ultra-high-resolution images of microvessels in a subject administered a microbubble contrast agent, as described in, for example, U.S. Patent Application Publication No. 16 / 617,628, and the contents of said U.S. Patent Application Publication are included in the contents of this application by reference. Generally, high resolution or ultra-high resolution refers to a resolution that is improved compared to the resolution achievable by the imaging system. For example, a high-resolution or ultra-high-resolution ultrasound image can refer to an image having a resolution higher than the diffraction limit.

[0015] The method includes the step of inputting ultrasound data into a computer system, as shown in step 102. In step 104, tissue motion can be removed from the ultrasound data. In step 106, the microbubble signal can be separated. In step 108, denoising processing can be performed on the microbubble signal data. In some embodiments, the location of the microbubbles can be identified in step 110, and the microbubbles can be tracked in step 112, which can be done, for example, after the velocity estimation method shown in Figure 11F below. In step 114, a microvascular image can be generated, processed in step 116, and displayed or stored in step 118.

[0016] In some embodiments, the step of inputting ultrasonic data into a computer system may include retrieving previously acquired ultrasonic data from memory or other data storage devices, which may be part of the computer system or separate from the computer system. In other embodiments, the step of inputting ultrasonic data may include acquiring such data using an ultrasonic system and inputting the acquired data into a computer system, which may be part of the ultrasonic system or separate from the ultrasonic system.

[0017] Ultrasonic data can be ultrasonic radio frequency ("RF") data or ultrasonic in-phase orthogonal ("IQ") data, etc. Generally, ultrasonic data has one or more spatial dimensions, which may include lateral dimensions, axial dimensions, elevation dimensions, and combinations thereof. For example, ultrasonic data may have two spatial dimensions, such as lateral and axial dimensions. Ultrasonic data may also have a temporal dimension, such as a slow time dimension (i.e., the temporal direction in which multiple ultrasonic signals are collected).

[0018] As described above, ultrasound data is acquired from subjects administered with a microbubble contrast agent. In some embodiments, multiple different microbubbles with different resonant ultrasound frequencies (e.g., microbubbles of different sizes) are used for imaging, and by selecting a specific ultrasound frequency (e.g., a frequency specific to either transmission or reception), only a subgroup of selected microbubbles can be imaged, resulting in ultrasound data containing a separated microbubble source. As another example, an ultrasound pulse with sufficient energy to rupture a specific number of microbubbles can be used, in which case the ruptured microbubbles release free gas bubbles from their microcapsules, generating ultrasound signals with different amplitudes than the intact microbubbles. This effectively creates a separated microbubble source that can be used for high-resolution imaging.

[0019] Microbubble signals can be obtained from both linear and nonlinear components of ultrasound. The linear component is typically the fundamental frequency of the applied ultrasound, while the nonlinear component can be the harmonic frequency, the fundamental frequency, or both. For example, the nonlinearity resulting from amplitude modulation imaging can be the fundamental frequency.

[0020] The ultrasound data is processed to remove tissue motion, as shown in step 104. During in vivo imaging, there may be transducer motion, tissue motion generated by the circulatory system (e.g., heartbeat motion, arterial pulsation), and tissue motion generated by the respiratory system. Since the amplitude of these motions is significantly larger than the size of the microvessels to be resolved, significant blurring occurs in the microvessel images, causing a decrease in the accuracy of blood flow velocity measurement. Therefore, the microbubble signals detected by ultrasound can be processed to remove the aforementioned tissue motion.

[0021] For example, ultrasound data can be processed to realign microbubble signals obtained from the same microvascular location but that are spatially misregistered. For instance, as shown in Figure 2, due to physiological movements such as respiration, tissue 202 moves to a different tissue location 204, and the target microvascular 206 moves to a different microvascular location 208. If this misregistration (i.e., positional displacement, movement) is not taken into account, it will cause errors in the final microvascular image.

[0022] To correct misregistration caused by tissue motion, image registration and tissue motion compensation methods can be used. For example, using image intensity-based image registration methods such as affine transforms, misregistered images can be returned to their original positions by translation, rotation, scaling, and shearing. Other usable image registration methods include image feature-based image registration methods, spectral-phase-based image registration methods, and transformation model-based image registration methods. As another example, using ultrasonic speckle tracking-based methods such as 2D normalized cross-correlation, motion vectors between a reference signal and the moving tissue and microvascular signals can be estimated. Subsequently, the signals can be realigned by reversing the motion based on the estimated motion vectors. As yet another example, a 2D phase shift map can be estimated between the reference image and the post-motion image, which can then be applied to the post-motion image in the Fourier domain, and the motion can be corrected by performing an inverse Fourier transform. As a non-limiting example, microbubble signal registration processing is performed between each frame of the ultrasonic microbubble signal and the reference frame of the microbubble signal (e.g., the first frame of the imaging sequence). Because ultrasound data is acquired at high frame rates, estimated tissue motion can generally be assumed to be smooth in the temporal direction, without sharp peaks that would suggest physiologically unrealistic, high-speed tissue motion. This characteristic of tissue motion can be used to suppress false peaks in tissue motion estimation and minimize misregistration.

[0023] Referring again to Figure 1, after tissue motion has been removed from the ultrasound data, the microbubble signal is separated from the ultrasound data as shown in step 106. Generally, the separation of the microbubble signal involves separating the microbubble signal from the background signal, such as the tissue signal and the signal from stationary microbubbles whose position does not change between acquisition frames (e.g., the signal when the microbubbles do not move between frames). In some embodiments, the microbubble signal can be separated using methods such as inter-frame signal difference, high-pass filtering in the time direction of the signal, or singular value decomposition ("SVD") based filtering.

[0024] As an example, time-based high-pass filtering can be used to isolate the microbubble signal. In such an embodiment, a cutoff frequency lower than the time frequency of the moving microbubble signal and higher than the time frequency of the stationary background signal can be used to filter the ultrasonic data and isolate the microbubble signal.

[0025] As another example, SVD-based filtering can be used, which may employ a singularity cutoff to separate background signals (e.g., tissue signals or non-moving microbubble signals, typically projected to lower singularities) from separated moving microbubble signals (typically projected to intermediate to higher singularities). As an example, to implement SVD-based filtering for extracting microbubble signals, the block-adaptive SVD filter described in the pending U.S. Patent Application Publication 2019 / 0053780 can be used. All the contents of that document are incorporated by reference into this application. As yet another example, to implement SVD-based filtering for extracting microbubble signals, the fast SVD filter described in the pending U.S. Patent Application Publication 2018 / 0220997 can be used. All the contents of that document are incorporated by reference into this application.

[0026] Optionally, the separated microbubble signals can then be denoised, as shown in step 108. The microbubble microvascular imaging technique disclosed herein is based on a power-based cross-correlation map of the microbubble signals and the PSF. In this way, denoising each frame of the microbubble signals can help reduce noise in the cross-correlation calculations, thereby improving the robustness of the microvascular imaging.

[0027] Generally, noise has characteristics similar to microbubble signals, and when noise is more intense and microbubble signals are less intense in deeper regions of tissue, it can become difficult to distinguish between the two. As a result, noise signals may be mistakenly labeled as microbubble signals. This can occur, leading to inaccurate depictions of microvessels and inaccurate velocity estimations.

[0028] As an example, noise reduction can be implemented using intensity-based thresholding. This method is relatively accurate when it can be assumed that the microbubble signal is more intense than the background noise signal. For example, a considerable amount of background noise can be suppressed by suppressing pixels with intensity values ​​lower than a selected value (e.g., values ​​from -30 dB to the maximum intensity value in the current field of view). However, this method may become inaccurate in areas where the microbubble signal resembles noise (e.g., deep tissue areas). Furthermore, the threshold must be carefully selected to avoid mistakenly excluding the microbubble signal or preserving an excessive amount of noise.

[0029] As another example, denoising of microbubble signals can be performed based at least partially on the spatiotemporal information contained in the microbubble signal. Since microbubbles move due to blood flow, the movement of microbubbles is a deterministic phenomenon that can be continuously tracked across multiple acquisition frames, whereas noise phenomena are random and do not show any track-like features when observed across multiple acquisition frames. The above-mentioned difference between microbubbles and noise can be used to achieve robust noise suppression in the spatiotemporal domain. As an example, a non-local mean ("NLM") denoising filter can be applied to the original microbubble data containing noise.

[0030] Another advantage of the spatiotemporal denoising filter described above is that, because the denoising process is performed in the spatiotemporal domain, little to no spatial blur occurs in the underlying microbubble signal. Other denoising methods (e.g., conventional Gaussian smoothing, Gaussian spectral apodization, wavelet thresholding, or iterative total variation ("TV") minimization) can be used in the spatiotemporal domain to achieve similar denoising effects. In some embodiments, axial-time microbubble signal data can be used for denoising, and in several other embodiments, lateral-time data or complete axial-lateral-time 3D data can also be used for denoising.

[0031] After denoising the microbubble signal data, a microvascular image can be generated based on the power-based cross-correlation technique of the present invention, as described in detail in this disclosure. The microvascular image can be displayed to the user as shown in step 118, or it can be stored for later use, for example, later analysis. In some embodiments, microvascular morphological measurements (e.g., vascular density and vascular torsion) can be estimated from the microvascular image. As another example, microvascular hemodynamic measurements (e.g., blood flow velocity and blood flow volume) can be estimated from the microvascular image. For example, a microvascular image can be presented by superimposing a B-mode image of a target tissue onto the microvascular image, or by placing such B-mode image side by side. Alternatively, a microvascular blood flow velocity image can be presented by superimposing a B-mode image of a target tissue onto a microvascular blood flow velocity image with color-encoded blood flow directions, or by placing such B-mode image side by side. In such embodiments, more than two blood flow directions can be indicated using multiple blood flow direction hues.

[0032] In step 110, conventional microbubble localization can be performed, and such microbubble localization may include determining the position of the microbubbles in each time frame of the microbubble signal data. For example, the center position of each separated microbubble signal can be determined so that the movement of the microbubbles can be tracked over time. Using this localized microbubble center position, ultra-high resolution microvascular images can be constructed, and the movement of the microbubbles can be tracked to calculate hemodynamic measurements such as blood flow velocity. Conventional ultra-high resolution imaging methods use the center position of the microbubbles It relies on specifics.

[0033] This disclosure provides a system and method for generating images with significantly improved spatiotemporal resolution by directly using a correlation map and using MB signals and the PSF to the power of N to improve the resolution of the cross-correlation map, without using the center position of the microbubbles. Below, with reference to Figure 3, a flowchart outlining each step of a non-limiting example of the cross-correlation-based method described above is shown. In step 310, ultrasound data can be accessed, for example, from an image archive, or ultrasound data can be acquired from a target tissue administered with microbubbles (MB) or another suitable ultrasound contrast agent. In step 320, the MB signal can be separated from the tissue signal by processing the ultrasound data using, for example, tissue clutter filtering techniques to detect the MB signal. In step 330, the signal-to-noise ratio (SNR) can be improved using preprocessing to prepare data for high-resolution imaging. In step 340, the PSF can be determined and an MB image can be generated. In step 350, the MB image and the PSF to the power of N can be determined and used to improve the resolution of the cross-correlation map, where N is the order of the exponentiation. The power of N of the MB image and PSF can determine the sharpness of the final microvascular image. In step 360, a cross-correlation map between each MB image and the system's PSF can be generated. In step 370, by integrating the power-based cross-correlation maps described above across a selected number of frames, a high-resolution, high-contrast image of the microvessels can be generated. Using the individual MB trajectories in the above spatiotemporal cross-correlation map, the estimated migration velocity and estimated migration direction of the MB can be provided, which can be used in step 380 to generate a hemodynamic image of blood flow.

[0034] The method shown in Figure 3 can be implemented using a block processing method (processing is performed block by block). This approach allows for separate microvascular image processing of smaller data subsets, and by combining all of these subsets, improved performance can be generated for microvascular images.

[0035] Continuing to refer to Figure 3, in step 310, ultrasound data can be accessed or acquired from a region of interest (ROI) of the target tissue to which an ultrasound contrast agent has been administered. Any suitable contrast agent can be used, including but not limited to MB, phase change droplets, nanodroplets, gas vesicles, etc. The ultrasound data can be in any suitable format, including but not limited to ultrasound RF data, ultrasound IQ data, ultrasound envelope data, ultrasound B-mode data, etc. The ultrasound data can be acquired using any suitable detection sequence, including but not limited to focused beamline imaging (line-by-line imaging), wide-beam multiline acquisition imaging, composite plane wave imaging, composite divergent wave imaging, synthetic aperture imaging, etc. The ultrasound data can be acquired based on linear ultrasound components, such as basic imaging, or based on nonlinear ultrasound components, such as pulse inversion ("PI") imaging, amplitude modulation ("AM") imaging, pulse inversion amplitude modulation ("PIAM") imaging, or a combination thereof. In general, acquired ultrasound data can have one or more spatial dimensions, which may include lateral, axial, and elevation dimensions, and combinations thereof. Ultrasound data may also have temporal dimensions, such as a slow time dimension (i.e., the temporal direction for acquiring multiple ultrasound frames or volumes).

[0036] In step 320 of Figure 3, MB signal detection and separation of the MB signal from the tissue signal are performed by extracting the signal representing blood flow information and excluding the background signal representing tissue clutter and noise. This can be done. Blood flow imaging can be enhanced by the signal separation described above. Generally, MB flows in the blood vessels together with blood cells and differs from tissue signals in terms of spatiotemporal characteristics, so MB signals can be extracted by applying appropriate clutter filtering techniques for blood flow imaging. In several non-limiting examples, MB can be extracted using interframe signal difference methods, high-pass tissue clutter filtering, SVD-based tissue clutter filtering, regression-based tissue clutter filtering, Eigen-based filtering, etc. A signal can be obtained. For example, in some embodiments such as CEUS, which can detect a signal from the MB using the nonlinear response of the MB and suppress tissue clutter, clutter filtering may not be necessary.

[0037] A high signal-to-noise ratio (SNR) of the obtained microvascular signal (MB) allows for the generation of optimized microvascular images. A low SNR can lead to false detection of MBs, potentially resulting in artifacts in the final microvascular image. The extracted MB data can be preprocessed in step 330 of Figure 3, which may include denoising, enhancement, and equalization of the MB signal. In one non-limiting example, an intensity threshold can be set to remove signals below a set intensity threshold as noise from the MB data, thereby providing enhanced MB data with a suppressed noise background. A fixed or adaptive threshold can be used. In another non-limiting example, the SNR of the MB data can be improved by denoising based on the spatiotemporal characteristics of the MB data. Typically, the movement of MBs across multiple ultrasound frames is constant, and such movement differs from the behavior of random noise; therefore, any denoising filter that operates on spatiotemporal data can be applied. Such noise reduction filters may include, but are not limited to, non-local mean (NLM) noise reduction filters, spatiotemporal low-pass filters of any form, smoothing filters, median filters, Savitsky-Goley filters, non-local mean filters, amplitude threshold discrimination methods, and the like. Since ultrasonic signal intensity can be spatially dependent due to the effects of time gain compensation (TGC), tissue attenuation, and beamforming processes, intensity equalization processing can also be used to equalize the spatially dependent intensity of the MB signal for subsequent processing. Non-limiting examples of such techniques, such as the use of noise distributions to equalize the MB signal, are described in U.S. Patent Application Publication No. 2019 / 0053780, the contents of which are incorporated herein by reference.

[0038] In another non-limiting example, the noise floor of an ultrasound image can be measured from an ultrasound system and used as a spatially variable threshold to suppress noise in the MB image. Background noise is typically spatially variable, and the noise level increases in deeper regions as the TGC used increases. In some configurations, effective noise suppression can be provided using a spatially variable threshold corresponding to the spatially variable noise level. By using the noise floor of the ultrasound system as a spatially variable threshold, pixels in the MB image that fall below that threshold can be considered noise and removed from the image. The spatially variable threshold can be controlled by applying a scale to the noise floor image, as shown below.

[0039]

number

[0040] Here, MB denoised (x,y) is the MB image after denoising, original (x,y) is the original MB image, and N(x,y) is the spatially fluctuating noise floor. α is a scaling factor that determines the overall threshold level, and α can be any positive value. If the background noise follows a normal distribution and the noise floor is estimated as the standard deviation of the noise, setting α=1 will exclude 68.2% of the noise pixels, setting α=2 will exclude 95.4% of the noise pixels, and setting α=3 will exclude 99.8% of the noise pixels. In practice, a larger scaling factor can suppress more noise better and increase the reliability of the saved MB image, but there is a trade-off in that low-intensity MB images will be removed.

[0041] Referring to Figures 12A to 12D, a non-limiting example of spatially dependent threshold discrimination based on noise floor estimation is shown. This noise floor can be measured or estimated using methods already disclosed in U.S. Patent Application Publication No. 2019 / 0053780, and all the contents of that document are included in the contents of this application by reference.

[0042] In one non-limiting example, noise data can be collected by turning off ultrasound transmission (e.g., setting the acoustic output to zero) and performing reception only using the same system configuration and imaging sequence as that used for actual blood flow signal acquisition (including identical transducer, TGC gain settings, receiver filter, beamforming settings, etc.).

[0043] In a non-limiting example, if there are ultrasonic systems where turning off the acoustic output is impossible or undesirable, noise data can be collected by transmitting the sound into the air at a minimum acoustic output. High-intensity echoes reflected from the air can be removed by clutter filtering, allowing for the extraction of pure noise data. Any suitable clutter filter can be used to remove high-intensity echoes from the air, such as high-pass tissue clutter filtering, SVD-based tissue clutter filtering, regression-based tissue clutter filtering, or eigenvalue-based filtering. From the collected noise data, the noise floor can be estimated as the variance of noise measurements in the time domain, the noise standard deviation, the noise mean, and other statistical properties.

[0044] Referring to Figures 12E and 12F, a non-limiting example of residual noise suppression after noise reduction using one of the above denoising methods is shown, where residual noise may still remain in the MB image. Assuming that residual noise is typically randomly distributed within the MB image and has a small spatial size, while the true MB signal is typically more consistent and relatively larger in size, a thresholding process based on the size of isolated objects in each MB image can be applied to remove objects smaller than a specific size threshold. A non-limiting example of MB images before (Figure 12E) and after (Figure 12F) object size thresholding is shown, and both figures demonstrate that the residual noise is better suppressed. To optimally suppress background noise in microvascular imaging, one of the above denoising methods can be applied individually, or a combination of several can be applied. Noise suppression methods can be applied to any conventional localization and tracking-based ultra-high-resolution technique, thereby achieving better denoising.

[0045] Subsequently, as shown in step 360 in Figure 3, the correlation or convolution between each frame of MB data and the obtained point spread function (PSF) can be calculated to generate a cross-correlation coefficient map. The PSF obtained in step 340 of Figure 3 may depend on the imaging settings, ultrasound system, center frequency, ultrasound pulse length, etc., and can be generated in any appropriate manner. This is possible. In some configurations, as described below, multiple different time frames of the cross-correlation map can be used in the velocity estimation method.

[0046] Determining the PSF may involve retrieving a pre-calculated or estimated PSF from memory or other data storage devices, or calculating or estimating the PSF as needed.

[0047] In a non-restrictive example, PSF can be obtained from simulation. PSF can be simulated, for example, based on a multivariate Gaussian distribution or other possible theoretical models. Such distributions are useful because they typically allow for ultrasonic lateral resolution coarser than axial resolution, although it is also possible to simulate PSF using other distributions.

[0048] As another non-limiting example, it is also possible to obtain PSF from measurements of very small point-like objects, such as objects significantly smaller than the wavelength of ultrasound. For example, the signal measured from microbubbles can be used as an approximation of the PSF of an ultrasonic system. PSF can be measured from point targets using the same ultrasonic system and imaging settings.

[0049] In another non-limiting example, PSF can also be obtained from individual isolated MBs derived from the same MB dataset for microvascular imaging. By selecting and combining multiple individual MB signals, a smoother and more accurate PSF can be generated. The PSF can be in any appropriate data format corresponding to the MB data used, such as RF data, ultrasound IQ data, ultrasound envelope data, ultrasound B-mode data, etc. A cross-correlation map can be generated by applying normalized or unnormalized correlation calculations. Correlation or convolution can be performed on a spatially interpolated version of the MB data, thereby obtaining a finer sample size. For each MB frame, true MB signals may show high correlation coefficients, while noise backgrounds or less relevant signals may have low correlation coefficients in the cross-correlation map. Thus, the cross-correlation map can provide information about the presence and location of MBs.

[0050] To improve the resolution of the cross-correlation map, the MB data to the power of N and the PSF to the power of N can be determined in step 350 of Figure 3 before calculating the cross-correlation map. Here, N can be any positive value, and in several non-limiting examples, N can be 1, 2, 3, 4, etc. When N=1, the power of the data is the same as the original data. One advantage of performing correlation processing on the MB data to the power of N and the corresponding PSF to the power of N is that a sharper MB target can be obtained for each MB in the cross-correlation map. In some configurations, the powers of the MB data and the powers of the PSF can be equal or different. In several non-limiting examples of the present disclosure, these two powers are set to the same value in the examples below.

[0051] Referring to Figure 4A, the resolution improvement of the N-th power method can be explained using an unrestricted example one-dimensional model. This one-dimensional model assumes that the original MB signal 410 is a one-dimensional Gaussian distribution, and therefore the PSF is also a one-dimensional Gaussian distribution. The cross-correlation or convolution of Gaussian distributions (MB signal and PSF) is also a Gaussian distribution, but with a larger width. Applying correlation calculations directly to the original MB data generates a cross-correlation map with worse resolution than the original MB signal, for example, the first-order correlation 420. In contrast, taking the N-th power cross-correlation of the Gaussian distribution can provide a sharper distribution curve with a smaller full width at half maximum (FWHM). As shown in the unrestricted example graph in the same figure, the 4th-order correlation 430, the 6th-order correlation 440, the 10th-order correlation 450, and the 20th-order correlation As can be seen from 460, the sharpness of the cross-correlation distribution curve gradually increases as the power degree N increases.

[0052] Referring to Figure 4B, an example of an FWHM graph is shown, which illustrates the quantitative relationship between the FWHM of the cross-correlation distribution and the power order N. The original line 510 shows the FWHM of the original MB distribution, while the cross-correlation line 520 shows the FWHM of the cross-correlation distributions obtained by various powers of N. When N > 2, the cross-correlation distribution begins to become sharper than the original MB distribution. In the Gaussian model, the FWHM can be improved at a rate of 1 / √N. The method disclosed herein can be applied to data of two dimensions or three or more dimensions.

[0053] Referring to Figure 5, in some configurations, the power of the MB data and the power of the PSF can be different. Figure 5 shows several non-restrictive examples of resolution improvement with power-based methods where the power of the MB data and PSF are different. As shown in the graphs of the non-restrictive examples in the same figure, the cross-correlation distribution curve gradually becomes sharper with increasing power of the MB data and the PSF, without requiring them to be the same.

[0054] Referring to Figures 6A-6D, non-restrictive example images obtained using the N-th power method at different power orders are shown. Figure 6A shows the original 2D MB signal. The cross-correlation map obtained by raising it to the first power is shown in Figure 6B, the cross-correlation map obtained by raising it to the fourth power is shown in Figure 6C, and the cross-correlation map obtained by raising the original MB signal to the eighth power is shown in Figure 6D. As the power order increases, the target distribution of each MB becomes smaller, and consequently, the sharpness of the resulting cross-correlation map can also be increased. In conventional ultra-high-resolution imaging, only the center point of each MB is labeled, resulting in a very sparse image containing many non-zero pixels, as shown in the example in Figure 6A. Conventionally, obtaining a visually smooth ultra-high-resolution image sometimes required integrating millions of MBs, which can necessitate very long acquisition times. Such improvements in resolution using conventional methods could not be achieved without very long acquisition times.

[0055] In contrast, the power-law correlation method disclosed herein allows the resolution to be changed by selecting the power. Each sharpened MB in the power-law correlation map can still have a large number of high-resolution pixels. The number of MBs required in the integration step to form a visually smooth, high-resolution microvascular image can be reduced. Data acquisition time can be shortened, and high spatial resolution can be achieved without requiring long acquisition times.

[0056] Referring to Figures 7A-7D, threshold-sharpened images of a non-restrictive example of cross-correlation are shown. In some configurations, to further enhance image sharpness, threshold discrimination can be applied to the resulting cross-correlation map to exclude correlation coefficient values ​​below a selected threshold. Using the thresholding method, small targets can be generated in the cross-correlation map for each individual MB. Figures 7A-7D show several non-restrictive examples of threshold discrimination on cross-correlation maps, and in these figures, the larger the threshold, the smaller the generated MB targets become. Therefore, the threshold can control the final resolution of the microvascular images.

[0057] In some configurations, denoising can be applied to remove suspicious targets in the cross-correlation map. Since genuine MB signals typically have a well-developed Gaussian shape, and noise and less relevant signals appear as irregular morphologies in the cross-correlation map, the well-developed Gaussian shape of the genuine MB signals can be utilized, potentially reducing the target size. As described above, these unwanted signals can be identified and removed from the correlation map based on their target morphology or target size. In a non-limiting example, thresholding can be applied to remove signals from a specific threshold. By removing smaller objects and storing larger objects in a cross-correlation map, the reliability and certainty of MB detection can be improved.

[0058] Referring to Figure 8, a non-limiting example of generating high-resolution, high-contrast microvascular images in step 370 of Figure 3 is shown, which is achieved by integrating a cross-correlation map over a selected number of frames. Conventional power Doppler blood flow images 820 are generated by directly integrating the original MB signal over a selected number of frames 810, and their resolution may be limited by the ultrasound wavelength and imaging settings. In a non-limiting example, the above integration can be performed based on a power-based cross-correlation map 830, which can generate microvascular images 840 with resolution beyond the limits of the imaging system. The resolution of the final microvascular image can be controlled by the power order and threshold discrimination processing of the cross-correlation map. In some configurations, a user interface can be provided to the user to allow the user to select or adjust the final resolution to meet the requirements of clinical applications.

[0059] In in vivo applications, tissue movement or ultrasound probe movement during scanning can cause problems. This can result in blurring of the final microvascular images, degrading the performance of this method. In some configurations, motion registration processing can be introduced to remove tissue motion and avoid blurring. Image registration processing can be performed based on motion estimation from the original acquired ultrasound data and / or clutter-filtered MB data. Any suitable image registration algorithm can be applied, including, but not limited to, global or local cross-correlation methods, global or local phase-correlation-based methods, and global or local optical flow methods.

[0060] Referring to Figure 9, multiple non-limiting images of two intracellular microvascular images obtained by the method disclosed in this application using a data length of approximately 2 seconds are shown in comparison with a conventional contrast-enhanced blood flow image. It is shown that the pancreatic microvascular image 920 is improved compared to the conventional pancreatic blood flow image 910, and that the liver microvascular image 940 is also improved compared to the conventional liver blood flow image 930.

[0061] In some configurations, additional image quality control can be applied before integration to suppress false MB signals and noise. True MB targets constantly move across multiple adjacent frames, which can result in continuous MB trajectories in the 3D spatiotemporal matrix of the correlation map (e.g., the matrix shown in image 830). In contrast, random noise or unrelated signals do not move continuously and may appear in the spatiotemporal matrix as isolated targets rather than long trajectories. The longer the trajectory, the more reliable the signal is as a valid MB. Therefore, thresholding can be performed in the 3D spatiotemporal matrix to exclude targets smaller than a selected threshold in size or length, thereby improving the reliability and certainty of MB detection in microvascular imaging.

[0062] Referring to Figure 10, several non-limiting examples of images with image quality control adjusted based on target size threshold discrimination are shown. Improved noise suppression for obtaining high-resolution images can be observed as the removal of objects smaller than the selected number of pixels. Image quality control can also be performed on the integrated 2D microvascular image obtained after the above-mentioned removal. In one non-limiting example, the integrated 2D image can be generated from a specific number of correlation maps, such as a correlation map of 10 frames, and threshold discrimination processing can be applied to exclude small targets, preserving MB trajectories in this 2D microvascular image. In another non-limiting example, the same process can be repeated, for example, by applying the same process to the correlation map of the next 10 frames, with or without overlap. By combining or integrating all of these separately threshold-discriminated microvascular images, the final microvascular image can be obtained.

[0063] Referring to Figure 11, (A) shows a non-restrictive example of MB trajectories in three-dimensional spacetime, which has two spatial dimensions and one temporal dimension. Each isolated MB trajectory in this three-dimensional spacetime is the spatial movement trajectory of the MB over time. In some configurations, MB movement trajectory analysis can be used to measure the flow velocity and direction of movement of each MB based on the three-dimensional spacetime matrix, as shown in step 380 of Figure 3. For example, high-resolution images of hemodynamic characteristics of blood flow (e.g., hemodynamic images) such as flow velocity, flow direction, and flow vector can be generated.

[0064] In some configurations, before velocity estimation, image quality control processing can be applied to the 3D spacetime to remove short or small trajectories and retain only long or large trajectories, thereby improving the robustness of velocity estimation. Figure 11(B) schematically shows an example MB trajectory in 3D spacetime. In Figure 11(B), the MB moves from the first frame 1115 to the Nth frame 1125, and each ellipse 1110 represents the MB in each of the multiple frames of the cross-correlation map. These MB ellipses 1110 constitute a single continuous trajectory in spacetime. The orientation and length of the trajectory can be determined by the MB's moving speed and direction of movement. By measuring the orientation and length of the trajectory, the MB flow velocity can be accurately estimated. Any suitable method can be used to measure the orientation and length of the trajectory.

[0065] In a non-restrictive example, the MB movement velocity can be measured for each MB trajectory by applying a fit to all discrete samples / pixels on that trajectory in 3D spacetime. Any fit can be applied, including, but not limited to, linear fit, linear regression, spline fit, square fit, hyperbolic fit, etc. Figure 11(C) shows an example of MB trajectory samples / pixels in 3D spacetime, and a linear fit to all of these trajectory samples can provide quantitative measurement results (including magnitude and direction) of the velocity vector. In a non-restrictive example, the slope of the fit line in the time direction can provide the estimated magnitude of the MB velocity, and the angle of the fit line in the spatial domain provides information about the MB movement direction.

[0066] In some configurations, MB trajectory samples can be projected onto a two-dimensional axial time plane as shown in Figure 11(D), or onto a two-dimensional lateral time plane as shown in Figure 11(E). These projections can be performed separately, and fitting the samples in both two-dimensional planes provides separate velocity measurements along the axial and lateral directions. The fitted lines in Figures 11(D) and (E) are fitted lines based on MB trajectory samples. The slopes of these fitted lines give the estimated MB velocity in both directions. One non-limiting example for further improving the robustness of velocity estimation is to apply weighted fitting to discrete samples / pixels of the MB trajectory. Any suitable weighted fitting can be used, including, but not limited to, weighted linear fitting, weighted linear regression, weighted spline fitting, weighted squared hyperbolic fitting, etc. Any suitable weights can be applied to the trajectory samples. In a non-restrictive example, trajectory samples can be weighted by their own cross-correlation coefficient, or by the k-th power of that cross-correlation coefficient, where k can be any value.

[0067] For each MB trajectory, the fitted correlation coefficient (R) and the coefficient of determination (R) are calculated. 2), the average of the fitting error The mean, standard deviation, or variance, or other statistical properties, can be used as indicators of fit robustness. Further image quality control for velocity estimation can be applied by removing MB trajectories with low fit robustness. Low fit robustness is associated with a low correlation coefficient (R) or coefficient of determination (R). 2 ) or suggested by a large fitting error, etc. To identify the desired fit robustness, set a threshold and identify MB trajectories that do not meet the set threshold and have low fit robustness, MB trajectories with a low correlation coefficient (R), and coefficient of determination (R 2 This may include removing MB trajectories with low accuracy or large fitting errors. In some configurations, velocity information (e.g., magnitude or direction of velocity) measured for each trajectory (each corresponding to one MB) can be stored spatiotemporally, and this can then be used to combine the velocity values ​​of all MB trajectories passing through the same microvessel to generate a final microvessel velocity image.

[0068] In some configurations, the measured velocity information of each trajectory can be assigned to a spatiotemporally shrinked trajectory, thereby further improving the spatial resolution of the final microvascular velocity image. In some configurations, the final microvascular velocity image can be generated by taking a simple or weighted average of all velocity trajectories in the spatiotemporal direction over time. For weighted averaging, velocity trajectories can be weighted using any possible weighting method, including the cross-correlation coefficient of the trajectory samples themselves, the fitted correlation coefficient (R), and the coefficient of determination (R). 2), including, but not limited to, the mean, standard deviation, or variance of the fitting error. In several non-limiting examples, microvascular morphological measurements (e.g., vascular density and vascular torsion) and microvascular hemodynamic measurements (e.g., blood flow velocity and blood flow volume) can be estimated from microvascular images. Microvascular images can be presented together with B-mode images of the target tissue, either superimposed or transparently superimposed, or with the B-mode images placed side by side. In another non-limiting example, microvascular velocity images, with or without color-encoded blood flow directions, can be superimposed or transparently superimposed on microvascular morphological images or B-mode images of the target tissue, or with the microvascular morphological images or B-mode images placed side by side, and with the microvascular velocity images presented together. More than two blood flow directions can be indicated using multiple blood flow direction hues.

[0069] In another non-restrictive example, each MB trajectory in 3D spacetime can be projected onto a 2D spatial domain, and the MB's velocity and direction of movement can be estimated using the projected area of ​​each MB trajectory. The faster the MB's velocity, the larger the projected area related to that MB; conversely, the slower the MB's velocity, the smaller the projected area may be. By calculating the projected area of ​​each separated MB trajectory, the estimated velocity of the MB can be obtained. In some configurations, the velocity can be estimated based on the relationship between the projected area of ​​the trajectory (A1) and the area of ​​the individual MB (area A0, for example, the area of ​​one ellipse 1110 in Figure 11(B)). The estimated velocity of the MB can be obtained from the ratio of A1 / A0.

[0070] In a non-limiting example, a fit 1120 can be applied to the feature points of each MB ellipse 1110, and the length and angle of the fit 1120 can be determined for velocity estimation. Any fit can be applied, including but not limited to linear fit, spline fit, square fit, hyperbolic fit, etc. Any suitable feature points of the MB ellipse can be used in this disclosure, including but not limited to the mean center of the MB ellipse, the weighted mean of the MB ellipse, the location of the maximum cross-correlation in the MB ellipse, the foci of the MB ellipse, the edges of the MB ellipse, etc.

[0071] In some configurations, the slope of the fitted line in the time direction can give the estimated magnitude of the MB velocity, while the angle of the fitted line in the spatial domain provides information about the direction of MB movement. In some configurations, MB feature points are initially placed in a two-dimensional axial-time plane and a two-dimensional... The data can be projected separately onto the lateral-time plane, and the fitting to these two two-dimensional planes provides separate measurement results for velocity in the axial and lateral directions.

[0072] To further improve the robustness of velocity estimation, weighted fitting of feature points can be applied. Any appropriate weighted fitting can be used, including, but not limited to, weighted linear fitting, weighted linear regression, weighted spline fitting, weighted squared fitting, and hyperbolic fitting. Any appropriate weight can be applied to the feature points. For each MB trajectory, the fitted correlation coefficient (R) and the coefficient of determination (R) are calculated. 2 The mean, standard deviation, or variance of the fitting error, or other statistical characteristics, can be used as indicators of fitting robustness. By removing MB trajectories with low fitting robustness, further image quality control for velocity estimation can be applied. Low fitting robustness is associated with a low correlation coefficient (R) or coefficient of determination (R). 2) or suggested by a large fitting error. In one non-limiting example, the length and direction of the trajectory can be roughly estimated using the distance and angle between the feature points of the first MB ellipse and the last MB ellipse 1110. In this way, it is sufficient to calculate only the feature points of the first MB ellipse and the last MB ellipse 1110, thereby reducing the computational load.

[0073] The steps of a non-limiting example of a velocity estimation method are shown in the flowchart of Figure 11F. In step 1150, the spatiotemporal microbubble trajectories are accessed from, for example, an image storage archive, or obtained from image data. In step 1155, individual microbubble trajectories can be separated. In step 1160, a fitting line can be applied to the microbubble trajectories. In some configurations, this fitting line can be applied to multiple projection planes, such as the lateral or axial projection planes described above. In step 1165, the orientation and length of the selected trajectories can be determined. In step 1170, optional image quality control can be performed by removing short trajectories or trajectories that do not fit the fitting line well, thereby improving the robustness of velocity estimation. Based on the determined orientation and length of the selected trajectories, the microbubble flow velocity can be estimated. Subsequently, in steps 1180 and 1185, microvascular images and hemodynamic images can be generated.

[0074] In some configurations, all of the field of view (FOV) data can be used. However, in some configurations, it may be advantageous to spatially divide the FOV data into multiple subsets of data with or without spatial overlap, and apply the high-resolution microvascular imaging method disclosed herein separately to each subset. By combining all of the microvascular images separated from all of these data subsets, a final high-resolution microvascular image can be obtained. Any appropriate data separation method can be used, and any appropriate weighting function can be applied. In a non-limiting example, the original data can be spatially divided into multiple smaller blocks with or without spatial overlap. In a non-limiting example, the original data can be spatially divided into multiple smaller blocks with spatial overlap, such as the data separation method described in U.S. Provisional Patent Application No. 62 / 975,515, and each data subset can be spatially weighted using a weighting function. The contents of that document are included herein by reference.

[0075] The separation of data into multiple subsets can be performed on the original data before tissue clutter filtering, on the MB data after tissue clutter filtering, and / or a combination thereof. When using the original data, an MB signal detection process, such as tissue clutter filtering, can be applied separately to each data subset. By dividing the original data into multiple subsets, an adaptive preprocessing unit can be applied, which can be unique to each subset. In a non-limiting example, an adaptive intensity threshold can be used so that signals below a certain threshold are considered noise and excluded from each subset of the MB data. The threshold for each subset of data can be adaptively determined based on local statistics of the MB data. In some configurations, the above threshold can be adaptively estimated based on the histogram of each data subset. In some configurations, the above threshold can be adaptively determined based on the total energy of each subset of the MB data.

[0076] In some configurations, the PSF can vary spatially, such as when ultrasonic focusing or beamforming apertures cause different PSF values ​​at different spatial locations. Different PSFs can be used for different subsets of data to compute spatially distinct cross-correlation maps. Using a unique PSF for each subset of data can be advantageous for achieving improved resolution and optimization for each localized dataset. Spatially distinct PSFs can be obtained through model fitting estimation, measurement, direct derivation from the same MB dataset, or a combination of these methods.

[0077] While this disclosure describes a method for MB imaging, it will be apparent to those skilled in the art that this method can be applied to ultrasound imaging using any other type of ultrasound contrast agent, such as phase-change droplets or nanodroplets, or to non-contrast ultrasound imaging without contrast agents. Although the method disclosed here is applied to two-dimensional imaging, it can also be applied to three-dimensional (3D) or higher-dimensional imaging. In a non-limiting example, the method can be easily extended to three dimensions, and a three-dimensional PSF can be generated to perform cross-correlation calculations in three dimensions. The method proposed in this disclosure can be combined with MB signal separation disclosed in U.S. Provisional Patent Application No. 62 / 861,580, in which the original MB data is separated into multiple subsets of sparsely distributed MB concentrations, thereby improving imaging performance. The contents of that document are incorporated herein by reference. The multiplication-based cross-correlation method disclosed herein can be combined with ultra-high-resolution methods such as, for example, location-based and tracking-based ultra-high-resolution methods (U.S. Patent Application No. 16 / 617,628), thereby improving the performance (accuracy and precision) of MB location determination in cross-correlation maps.

[0078] Figure 13 shows an example of an ultrasonic system 1300 capable of carrying out the method disclosed herein. The ultrasonic system 1300 comprises a transducer array 1302 having a plurality of separately driven transducer elements 1304. The transducer array 1302 can include any suitable ultrasonic transducer array, including linear arrays, curved arrays, phased arrays, etc. The transducer array 1302 can also include 1D transducers, 1.5D transducers, 1.75D transducers, 2D transducers, 3D transducers, etc.

[0079] When energized by transmitter 1306, a specific transducer element 1304 generates a burst of ultrasonic energy. The ultrasonic energy (e.g., echo) reflected from the object or subject and returned to the transducer array 1302 is converted into an electrical signal (e.g., an echo signal) by each transducer element 1304 and can be separately applied to receiver 1308 by a set of switches 1310. Transmitter 1306, receiver 1308, and switches 1310 operate under the control of controller 1312, which may comprise one or more processors. As an example, controller 1312 may comprise a computer system.

[0080] Transmitter 1306 can be programmed to transmit unfocused or focused ultrasound. In some configurations, transmitter 1306 can also be programmed to transmit divergent waves, spherical waves, cylindrical waves, plane waves, or a combination thereof. Transmitter 1306 can also be programmed to transmit spatially or temporally encoded pulses.

[0081] The receiver 1308 can be programmed to implement a suitable detection sequence for the current imaging task. In some embodiments, this detection sequence may include one or more of line-by-line scanning, synthetic plane-wave imaging, synthetic aperture imaging, and synthetic divergent beam imaging.

[0082] In some configurations, the transmitter 1306 and receiver 1308 can be programmed to achieve a high frame rate. For example, a frame rate corresponding to an acquisition pulse repetition frequency (PRF) of at least 100 Hz can be achieved. In some configurations, the ultrasonic system 1300 can sample and store at least 100 ensembles of echo signals in the time domain.

[0083] Scanning can be performed by setting each switch 1310 to its transmit position during a single transmission process according to a selected imaging sequence, thereby turning on the transmitter 1306 to energize the transducer element 1304. Subsequently, the switch 1310 is set to the receive position, and the successive echo signals generated by the transducer element 1304 in response to one or more detected echoes are measured and applied to the receiver 1308. The separate echo signals from each transducer element 1304 can be combined in the receiver 1308 to form a single echo signal.

[0084] The echo signal is sent to the processing unit 1314 for processing, or an image generated from the echo signal is processed. This processing unit 1314 can be implemented using a hardware processor and memory. For example, the processing unit 1314 can generate an image using the method disclosed herein. The image generated from the echo signal by the processing unit 1314 can be displayed on the display system 1316.

[0085] Referring to Figure 14, a block diagram of an example of a computer system 1400 capable of carrying out the method disclosed herein is shown. The computer system 1400 generally comprises an input unit 1402, at least one hardware processor 1404, a memory 1406, and an output unit 1408. Thus, the computer system 1400 is generally embodied by the hardware processor 1404 and the memory 1406.

[0086] In some embodiments, the computer system 1400 may be a controller or processing unit for the ultrasonic system. The computer system 1400 may also be embodied in some examples by a workstation, notebook computer, tablet device, mobile device, multimedia device, network server, mainframe, one or more controllers, one or more microcontrollers, or any other general-purpose or purpose-specific computer.

[0087] The computer system 1400 can operate autonomously or semi-autonomously, or it can read executable software instructions from memory 1406 or a computer-readable medium (e.g., a hard drive, CD-ROM, flash memory, etc.), or it can receive instructions from a user via input unit 1402, or it can receive instructions via input unit 1402 from other sources logically connected to the computer or device, such as other networked computers or servers. Therefore, in some embodiments, the computer system 1400 may also be equipped with any suitable device for reading computer-readable storage media.

[0088] Generally, the computer system 1400 is programmed to implement the methods and algorithms disclosed herein, or is otherwise configured. For example, the computer system 1400 may be programmed to generate images using the methods disclosed herein.

[0089] The input unit 1402 can, as desired, take any suitable shape or form for the operation of the computer system 1400, including the ability to perform tasks, process data, or select, input, or otherwise specify parameters corresponding to the operation of the computer system 1400. In some aspects, the input unit 1402 can be configured to receive data, such as data acquired using an ultrasonic system. Such data can be processed as described above to generate an image. Furthermore, the input unit 1402 can also be configured to receive any other data or information that is considered useful for generating an image using the method described above.

[0090] In processing tasks for operating the computer system 1400, one or more hardware processors 1404 can be configured to perform any number of post-processing steps on data received using the input unit 1402.

[0091] The memory 1406 can store the software 1410 and data 1412, such as data acquired using an ultrasonic system, and can also be configured to store and retrieve processed information, instructions, and data processed by one or more hardware processors 1404. In some aspects, the software 1410 may include instructions for image generation according to embodiments of the present disclosure.

[0092] Furthermore, the output unit 1408 can take on any shape or form as desired, and can be configured to display images (e.g., microvascular images), and can also display other desired information.

[0093] In some embodiments, any suitable computer-readable medium can be used to store instructions for performing the functions and / or processes described herein. For example, in some embodiments, the computer-readable medium may be temporary or non-temporary. Non-temporary computer-readable mediums may include, for example, magnetic media (e.g., hard disks, floppy disks, etc.), optical media (e.g., compact discs, digital video discs, Blu-ray discs, etc.), semiconductor media (e.g., random access memory ("RAM"), flash memory, electrically programmable read-only memory ("EPROM"), electrically erasable programmable read-only memory ("EEPROM"), etc.), any suitable medium that is not transient or loses any form of performance during transmission, and / or any suitable tangible medium. As another example, temporary computer-readable mediums may include signals on a network, in wiring, in conductors, in optical fibers, on circuits, or on any suitable medium that is not transient or loses any form of performance during transmission, and / or any suitable intangible medium.

[0094] While this disclosure describes one or more preferred embodiments, many equivalent, alternative, modified, and improved embodiments are possible in addition to those explicitly described, and it is clear that these embodiments fall within the scope of the present invention.

Claims

1. A method for performing high spatiotemporal resolution imaging of microvessels using an ultrasound system, (a) Accessing the ultrasound data acquired using the ultrasound system from the region of interest of a subject in which a microbubble contrast agent was present when the ultrasound data was acquired, using a computer system; (b) A step of generating microbubble signal data by separating the microbubble signals in the ultrasonic data from other signals in the ultrasonic data using the computer system, (c) A step of separating the individual microbubble trajectories in the microbubble signal data, (d) A step of determining the trajectory parameters of the separated microbubble trajectories, (e) A step of estimating the microbubble flow velocity based on the trajectory parameters obtained from the separated microbubble trajectories, (f) A step of generating a microvascular image with high spatiotemporal resolution based at least partially on the estimated microbubble flow velocity, has A method characterized by the following:

2. The step of determining the trajectory parameters includes applying a fitting line to each of the microbubble trajectories. The method according to claim 1.

3. The aforementioned trajectory parameter includes at least one of direction or length. The method according to claim 2.

4. At least one of the orientation or length is determined based on the movement speed and direction of the microbubbles between multiple different frames of the cross-correlation map. The method according to claim 3.

5. The method according to claim 2, wherein the slope of the fitting line with respect to time provides the magnitude of the microbubble velocity.

6. The angle of the fitting line in the spatial domain provides the direction of movement of the microbubbles. The method according to claim 2.

7. The aforementioned fitting line is a weighted fit of discrete samples of the separated microbubble trajectories, The aforementioned weighted fitting includes at least one of the following: weighted linear fitting, weighted linear regression, weighted spline fitting, weighted square fitting, and weighted hyperbolic fitting. The weights in the weighted fitting are weights based on the cross-correlation coefficients of the separated microbubble trajectories, or weights based on powers of the cross-correlation coefficients. The method according to claim 2.

8. The aforementioned fitting lines are fitted to the characteristic points of each microbubble trajectory. The aforementioned feature point includes at least one of the following: a microbubble ellipse, the mean center of the microbubble ellipse, the weighted mean of the microbubble ellipse, the position of maximum cross-correlation in the microbubble ellipse, the focus of the microbubble ellipse, and the edge of the microbubble ellipse. The method according to claim 2.

9. The process further includes a step of determining an index of fitting robustness of the aforementioned fitting line, The aforementioned fit robustness index is the fit correlation coefficient (R) of the fitted line and the coefficient of determination (R) of the fitted line. ２ ) and at least one of the mean, standard deviation, or variance of the fitting error of the aforementioned fitting line, The method according to claim 2.

10. Before estimating the microbubble flow velocity, the method further includes the step of removing the separated microbubble trajectories by removing those microbubble trajectories that do not reach a threshold of at least one of the following indicators: the fitting correlation coefficient (R) of the fitted line, the coefficient of determination (R²) of the fitted line, and the mean, standard deviation, or variance of the fitting error of the fitted line. The method according to claim 9.

Images