Phase velocity imaging using an imaging system

Abstract

Described herein are systems and methods for phase velocity imaging using imaging systems such as ultrasound systems, optical imaging systems (e.g., optical coherence tomography systems), or magnetic resonance imaging systems. Generally, systems and methods are described for constructing phase velocity images (e.g., 2D images, 3D images) from propagating mechanical wave data. The systems and methods described in this disclosure operate in the frequency domain and can be implemented using a single frequency or a band of selected frequencies. When there are multiple mechanical wave sources within the field of view, directional filtering can be performed to separate mechanical waves propagating in different directions. The reconstruction described below can be performed on each of these directionally filtered components.

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of U.S. Provisional Patent Application No. 62 / 740,730, entitled “PHASE VELOCITY IMAGING USING AN ULTRASOUND IMAGING SYSTEM,” filed October 3, 2018, which is incorporated by reference in its entirety.

[0002] Statement regarding federally sponsored research This invention was made with Government support under DK092255 awarded by the National Institutes of Health. The Government has certain rights in this invention. [Background technology]

[0003] Soft tissue pathology can cause changes in mechanical properties. Many methods for measuring the mechanical properties of soft tissue have been developed over the past two decades. One class of methods, called shear wave elastography ("SWE"), includes a variety of different implementations, including supersonic shear imaging, comb-push ultrasonic shear elastography, among others. These methods use acoustic radiation force to generate shear waves and ultrasound techniques to measure shear wave motion.

[0004] These methods have been employed in imaging a variety of tissues, including liver, breast, thyroid, skeletal muscle, kidney, and prostate. Liver fibrosis staging has been the primary and successful application of SWE techniques. SWE has also been widely used for characterization of breast, thyroid, and prostate tumors. In SWE applications related to breast, thyroid, and liver cancer imaging, inclusions have typically been found to be stiffer than normal surrounding tissue. These methods have been used to characterize lesions as benign and malignant, and malignant lesions have typically been found to be stiffer than benign lesions. In clinically implemented versions of SWE, several assumptions are made about the medium. These include that the material is elastic, homogeneous, isotropic, and incompressible. For imaging inhomogeneities of finite size, resolution is an important factor to consider.

[0005] Among the different methods used to implement SWE measurements, there are also numerous techniques to reconstruct the mechanical properties, especially the shear wave velocity. Most methods are based on time domain data and measure the time of flight of the shear wave in a local sense. Different features of the shear wave motion are used to estimate the shear wave velocity, including the peaks of the displacement or velocity waveforms. Furthermore, cross-correlation of the time signal can be used and the size of the window used can affect the shear wave velocity result.

[0006] Another way to analyze the data to estimate the shear wave velocity is to use a frequency domain approach. Typically, frequency domain evaluation of the shear wave velocity has been performed with phase gradient or Fourier transform based methods. These methods have been used to estimate the phase velocity dispersion due to the viscoelasticity of the material and the geometry. Most of these measurement techniques, as implemented, use limited spatial ranges in both the lateral and axial directions. To fill the gap, efforts have been made to extend these methods to produce images of the phase velocity, or viscoelastic properties such as shear wave attenuation and the viscosity of the material. Summary of the Invention

[0007] The present disclosure addresses the aforementioned shortcomings by providing a method for generating a phase velocity image from mechanical wave data acquired with an ultrasound system, accessing mechanical wave data acquired with an ultrasound system and generating phase velocity data therefrom, and then generating a phase velocity image from the phase velocity data.

[0008] Another aspect of the present disclosure is to provide a method for controlling an ultrasound system to acquire mechanical wave data. Optimal acquisition parameters of the ultrasound system are determined by (a) acquiring mechanical wave data with the ultrasound system over a range of frequency values, (b) reconstructing a phase velocity image from the mechanical wave data, and (c) analyzing the phase velocity image to determine whether a phase velocity plateau is present in the phase velocity image. If the phase velocity plateau is not present in the phase velocity image, the acquisition parameters of the ultrasound system are adjusted and steps (a)-(c) are repeated using different frequency values ​​until the phase velocity plateau is present in the phase velocity image. The adjusted acquisition parameters are stored as optimal acquisition parameters. The optimal acquisition parameters are associated with an optimal bandwidth of the mechanical wave to improve accuracy of the phase velocity measurement.

[0009] The foregoing and other aspects and advantages of the present disclosure will become apparent from the following description. In the description, reference is made to the accompanying drawings, which form a part hereof, and in which a preferred embodiment is shown by way of example. This embodiment does not, however, necessarily represent the full scope of the invention, and therefore, reference is made to the claims and this specification in order to interpret the scope of the invention. [Brief description of the drawings]

[0010] [Figure 1] 4 is a flowchart illustrating steps of an exemplary method for generating a phase velocity image from wave data acquired with an imaging system. [Diagram 2]1 is a flowchart illustrating steps of an exemplary method for generating a phase velocity image from wave data acquired with an imaging system using a local phase velocity imaging ("LPVI") algorithm. [Diagram 3] 1 is a flowchart illustrating steps of an example of an LPVI method for generating a phase velocity image from wave data. [Figure 4] FIG. 1 is a block diagram of an example system capable of implementing the methods described in this disclosure. [Diagram 5] FIG. 5 is a block diagram of exemplary hardware components that may be implemented in the system of FIG. 4. [Figure 6] FIG. 1 is a block diagram of an example ultrasound imaging system that can be implemented to acquire wave data. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0011] Described herein are systems and methods for phase velocity imaging using ultrasound or other suitable imaging systems, such as optical imaging systems or magnetic resonance imaging ("MRI") systems. Generally, systems and methods for constructing phase velocity images (e.g., 2D images, 3D images) from propagating wave data are described. The wave data can be guided wave data, shear wave data, or other such mechanical wave data. The systems and methods described in this disclosure operate in the frequency domain and can be implemented using a single frequency or a band of selected frequencies. When there are multiple acoustic radiation force pushes in the field of view or multiple wave sources, directional filtering can be performed to separate waves propagating in different directions. The reconstruction described below can be performed on each of these directionally filtered components.

[0012] The systems and methods described in this disclosure enable the reconstruction of phase velocity images with sufficient resolution to accurately depict inclusions or other objects on the order of less than 5 mm in diameter, objects of that size that have traditionally been difficult to resolve with previous time-domain based methods.

[0013] In some aspects of the present disclosure, systems and methods for local phase velocity imaging ("LPVI") are described. LPVI implements local wave number imaging with the goal of creating an image of the phase velocity of soft tissue. For example, a local estimate of the wave number of the propagating guided Lamb waves can be utilized. The change in wave number as a function of spatial location can be analyzed with these systems and methods. For example, new wave number components resulting from abrupt wave changes at structural discontinuities can be identified in the frequency-wave number spectrum. LPVI can accurately reconstruct 2D wave velocity maps (e.g., phase velocity maps) and provide good contrast between inclusions and surrounding background tissues without artifacts.

[0014] Phase Velocity Imaging Referring now to FIG. 1, a flow chart is shown illustrating steps of an exemplary method for phase velocity imaging using an imaging system. m (m,n) are accessed by the computer system in step 102 and then optionally directionally filtered into q distinct components in step 104 as desired or necessary. For example, if the wave data is acquired using an acquisition scheme in which multiple wave sources are directed in a single acquisition (e.g., by using multiple distinct push beams or other wave sources), directional filtering can advantageously reduce wave interference due to reflection and refraction. Additionally or alternatively, wave number filtering may also be applied to the wave data, as shown in optional process block 106.

[0015] Accessing the wave data can include retrieving the wave data from a memory or other data storage device or medium. In other examples, accessing the wave data can include acquiring such data with an imaging system and communicating or otherwise transferring the wave data to a computer system that may be part of the imaging system. The imaging system can be an ultrasound system, an optical imaging system (e.g., an optical coherence tomography system), an MRI system (e.g., by using magnetic resonance elastography techniques), or any other imaging system suitable for imaging propagating mechanical waves (e.g., shear waves, guided waves, other mechanical waves).

[0016] The wave data may be mechanical wave data, such as shear wave data, induced wave data, etc., and may indicate particle displacement, velocity, acceleration, or a combination thereof. The wave data may be 3D wave data (e.g., two spatial dimensions, one time dimension) or 4D wave data (e.g., three spatial dimensions, one time dimension). Although the method is described below with respect to processing 3D wave data (e.g., two spatial dimensions and one time dimension), it will be understood by those skilled in the art that the method can be easily extended to receive 4D wave data (e.g., three spatial dimensions, one time dimension) as input.

[0017] Different phase velocity estimation algorithms can then be employed in step 108 to calculate or otherwise generate phase velocity data, including lateral, axial, and elevation phase velocity components, i.e., v X,DF,q (m,n,f),v Z,DF,q (m,n,f), and v Z,DF,q(m,n,f). This step may isolate one frequency or a set of selected frequencies (e.g., a frequency band). In general, increasing the selected frequency may improve the quality of the reconstructed phase velocity image in the sense that the phase velocity can be made more uniform and the shape of the inclusions can be more accurately represented. Furthermore, the magnitude of the wave velocity (e.g., shear wave velocity) increases with increasing frequency, which allows better contrast to be obtained. The spatial resolution of the reconstructed phase velocity image is determined by the shorter spatial wavelength at higher frequencies. This is advantageous for resolving small inclusions with high contrast. In some cases, using a single frequency may not provide a proper signal-to-noise ratio. Using a frequency band may allow the overall spatial wavelength to be maintained, but the signal-to-noise ratio to be improved.

[0018] In some implementations, a moving patch of p pixels may be used to calculate the phase velocity. The phase velocity estimation algorithm may be implemented in a hardware processor and memory and may include a phase gradient ("PG") algorithm, a two-dimensional Fourier transform algorithm, a two-point continuous wavelet ("2P-CWT") algorithm, a multiple signal classification ("MUSIC") algorithm, an LPVI algorithm, or other suitable phase velocity estimation algorithm.

[0019] After the individual phase velocity components are estimated at a representative patch p, a window can be used for averaging in step 110. The averaging can be Anderssen-Hegland averaging using the methods described in R S Anderssen and M. Hegland, "For numerical differentiation, dimensionality can be a blessing!", Mathematics of Computation, 1999;68(227):1121-1141, P. Song et al., "Fast shear compounding using robust 2-D shear wave speed calculation and multi-directional filtering", Ultrasound Med.Biol., 2014;40:1343-1355, and U.S. Patent No. 9,622,711, which are incorporated herein by reference in their entireties.

[0020] Using this approach for a 2D-processing window, the final wave velocity V at the central pixel V(m,n) is given by V X and V Z and then using these maps in equations (1) and (2).

number

number

[0021] where h′ and r(i,j) are defined as follows:

number

[0022] and m is the lateral dimension of a pixel, n is the axial dimension of a pixel, and CC X and C.C. Zare the normalized cross-correlation coefficients along the x and z directions.

[0023] Phase velocity component (V X,AH,DF,q (m,n,f), V Z,AH,DF、q (m,n,f), V Y,AH,DF,q (m,n,f) can be combined in either two or three dimensions, as shown in step 112. As an example, the phase velocity components can be combined with the following equation:

number

number

[0024] If the wave data is directionally filtered, the directionally filtered components can be combined in step 114 to generate a phase velocity image. As an example, the directionally filtered components can be combined with weights b q They can be combined using a weighted sum of (m,n).

number

[0025] The above described method provides a phase velocity image that can be used in a variety of methods and applications. One advantage of a phase velocity image is that it can be reconstructed at known frequencies or using known frequency bands, which allows for standardization.

[0026] For viscoelastic materials, the phase velocity varies with frequency or dispersion. The phase velocity dispersion can be fitted to rheological models such as Kelvin-Voigt, Maxwell, Generalized Maxwell, Zener, or fractional models to estimate the viscoelastic parameters and storage and loss moduli. As a result, material property maps (e.g., viscoelastic parameter maps, mechanical parameter maps) can be optionally generated, as shown in step 314. These parameters provide an additional way to characterize soft tissues and other materials.

[0027] The reconstruction of phase velocity images described in this disclosure can be optimized using different metrics such as coefficient of variation (i.e., ratio of the mean to the standard deviation within a region), standard deviation or variance within a region, contrast between two different regions, contrast-to-noise ratio ("CNR") between two different regions, and resolution. Processing parameters that can be varied in an iterative or guided manner include patch size p, window size w, frequency or frequency band, or combinations of these parameters. The optimization process can be performed in an iterative or guided process to maximize a metric (e.g., contrast, CNR) or minimize a metric (e.g., coefficient of variation, standard deviation, variance, resolution). The size of the region or inclusion can be used to initialize the optimization process. Furthermore, the group velocity image and an examination of the center frequency of the shear wave can guide the selection of frequency and size processing parameters.

[0028] In combination with Anderssen-Hegland averaging, various one-dimensional methods for measuring phase velocities can be adapted to two and / or three dimensions to provide a measure of the true wave velocity in 2D or 3D, not just the wave velocity of lateral traveling waves. Furthermore, Anderssen-Hegland averaging can provide a level of robustness to the reconstruction.

[0029] In some implementations, the methods described in this disclosure can be implemented on different sets of wave data acquired when generating shear waves in different parts of the same volume and combining the resulting phase velocity images. For example, a first phase velocity map can be generated from wave data acquired when a push beam excitation is generated on one side (e.g., the left side) of the volume (e.g., using the left side of the ultrasound transducer to generate the push beam excitation). Then, a second phase velocity map can be generated from wave data acquired when a push beam excitation is generated on another side (e.g., the right side) of the volume (e.g., using the right side of the ultrasound transducer to generate the push beam excitation). These first and second phase velocity images can then be combined, such as by averaging the two images.

[0030] Phase gradient and phase velocity estimation As mentioned above, in some implementations, the phase velocity estimation algorithm can be a phase gradient ("PG") algorithm. In such methods, a phase gradient calculation is used to reconstruct the phase velocity of a local region defined as a patch p, pixel by pixel or spatial location by spatial location. To find the phase values ​​for the phase gradient calculation, a one-dimensional Fourier transform is performed in the time dimension of each recording of the wave motion at each location. The phase values ​​at each location and for each frequency or set of frequencies are extracted for later use in calculating the phase gradient. For each frequency of interest, the phase gradient

number

number

[0031] Two-dimensional Fourier transform (2DFT) As noted above, in some implementations, the phase velocity estimation algorithm may be a two-dimensional Fourier transform ("2DFT") algorithm. In such methods, the 2D Fourier transform of the spatiotemporal data v m (x,t) is calculated over a certain range containing p pixels. The spatial length is x T =px p Wherein, x p is the pixel width. Using a 2D Fourier transform provides a representation in the spatial and temporal frequency space V(k,f). For each frequency, a peak appears in the magnitude spectrum |V m (k,ω)|. The location of the peak is used to calculate the phase velocity using:

number

[0032] The location at which the phase velocity is calculated can also be found for each wavenumber k in a similar manner as was done for each frequency. These methods are further described in "Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes" by M. Bernal et al., J. Acoust. Soc. Am., 2011;129:1344-1354, which is incorporated herein by reference in its entirety.

[0033] When using short segments of only p pixels, zero padding may need to be performed. Additionally, wave propagation mode identification can be used by using a gradient-based methodology and fitting a polynomial curve to the data points to provide a more robust estimate of the dispersion curve, as described in P. Kijanka et al., "Robust phase velocity dispersion estimation of viscoelastic materials used for medical applications based on the Multiple Signal Classification Method," IEEE Trans Ultrason Ferroelectr Freq Control, 2018;65:423-439, which is incorporated herein by reference in its entirety.

[0034] Two-point continuous wavelet method (2P-CWT) As mentioned above, in some implementations, the phase velocity estimation algorithm can be a two-point continuous wavelet ("2P-CWT") algorithm. In such methods, a wavelet g(t), such as the Morlet wavelet, is composed of harmonics modulated by a Gaussian envelope centered around T0.

number

[0035] The mother wavelet g(t) can be expanded by α and delayed by τ as a daughter wavelet.

number

[0036] The wavelet transform is defined as follows:

number

[0037] where "*" denotes the complex conjugate. For a given separation of p spatial locations, the phase velocity can be estimated at a specified set of frequencies. p For two signals at different spatial locations separated by p pixels, we can compute the continuous wavelet, as follows, to calculate the correlation and the phase shift at a particular frequency from the unwrapped phase:

number

number

[0038] Further details of such methods are described in Q. Wu et al., "Measurement of interstation phase velocity by wavelet transformation," Earthquake Science, 2009;22:425-429.

[0039] Multiple signal classification (MUltiple SIgnal Classification: MUSIC) As mentioned above, in some implementations, the phase velocity estimation algorithm can be a multiple signal classification ("MUSIC") algorithm. In such methods, the phase velocity is calculated as a function of dimension x pThe phase velocity can be estimated for a segment of p pixels of . The phase velocity can be identified by fitting a polynomial curve to the data points using a gradient-based methodology to provide a more robust estimate of the dispersion curve. An exemplary MUSIC algorithm is described in Patent Application No. PCT / US2018 / 036050, which is incorporated by reference in its entirety, and in P. Kijanka et al., "Robust phase velocity dispersion estimation of viscoelastic materials used for medical applications based on the Multiple Signal Classification Method," IEEE Trans Ultrason Ferroelectr Freq Control, 2018;65:423-439.

[0040] Local Phase Velocity Imaging (LPVI) As mentioned above, in some implementations, the phase velocity estimation algorithm can be a local phase velocity imaging ("LPVI") algorithm. In such a case, the method described with respect to FIG. 1 can be modified as shown in FIG. 2, which shows a flow chart describing steps of an exemplary method for phase velocity imaging using an ultrasound system and an LPVI algorithm. The LPVI method can be used to obtain phase velocity images with or without wavenumber filtering. Anderssen-Hegland averaging can be used as an additional step to the LPVI step.

[0041] In general, the wave data v m(m,n) are accessed by the computer system in step 202 and then optionally directionally filtered into q distinct components in step 204 as desired or necessary. Additionally or alternatively, wave number filtering may also be applied to the wave data, as shown in optional process block 206. Accessing the wave data may include retrieving the wave data from a memory or other data storage device or medium. In other examples, accessing the wave data may include acquiring such data with an imaging system and communicating or otherwise transferring the wave data to a computer system that may be part of the imaging system. The imaging system may be an ultrasound system, an optical imaging system, an MRI system, or any other imaging system suitable for imaging propagating mechanical waves (e.g., shear waves, guided waves, other mechanical waves).

[0042] The wave data may be mechanical wave data, such as shear wave data, guided wave data, etc., and may indicate particle displacement, velocity, acceleration, or a combination thereof. The wave data may be 3D wave data (e.g., two spatial dimensions, one time dimension) or 4D wave data (e.g., three spatial dimensions, one time dimension). Although the method is described below with respect to processing 3D wave data, one skilled in the art will appreciate that the method can be readily extended to receive 4D wave data as input.

[0043] An LPVI algorithm is then employed in step 208 to calculate or otherwise generate phase velocity data, such as by calculating the phase velocities of the different directionally filtered components. For example, the directionally filtered components of the wave data can be input to the LPVI algorithm to generate an output as phase velocity data. This step may isolate one frequency or a selected set of frequencies (e.g., a frequency band). The directionally filtered components can then be combined in step 210 to generate a phase velocity image. For example, the directionally filtered components can be combined according to equation (5):

[0044] As mentioned above, LPVI algorithms operate in the "k-space" (frequency-wavenumber) domain where additional information about wave modes and wavenumber distributions can be obtained compared to the space-time domain. In general, the transformation of the space-time particle motion v(z;x;t) into the k-space domain is performed using a three-dimensional (3D) Fourier transform (F 3D ) is used to achieve this.

number

[0045] In the formula, V(k Z ,k X ,f) is the frequency and wave vector k Z and k X This is the 3D k-space representation obtained by

[0046] The frequency f corresponds to the time domain, whereas the wave number k Z and k X correspond to the spatial dimensions z and x, respectively. The obtained spectral data V(k Z ,k X , f), a wavenumber spectrum can be obtained for a particular frequency f0.

[0047] 3, a flow chart is shown illustrating steps of an exemplary method for generating phase velocity data using the LPVI algorithm. In LPVI, a transformation of the space-time particle motion v(z,x,t) can be performed step-by-step. The LPVI approach reconstructs local wave velocities (e.g., shear wave velocity, guided wave velocity) for a single frequency or a selected frequency band.

[0048] The method begins with transforming wave data into the frequency domain, as shown in step 302. In some cases, transforming wave data (e.g., spatiotemporal particle motion data v(z,x,t)) into the frequency domain can be achieved by a one-dimensional Fourier transform (F 1D ) can be used.

number

[0049] During the ceremony,

number

number

[0050] In other examples, the wave data is first transformed into the frequency-wavenumber domain (eg, using a 3D Fourier transform) and the resulting spectral data is directionally filtered, bandpass filtered, or both.

[0051] As an example of directional filtering, a directional filter can be applied to the 2D wavenumber domain to extract directional components, such as left-to-right ("LR") and right-to-left ("RL") propagating wave components, and remove wave interference. Axially propagating wave components can also be filtered out, k ZThis can be reduced by directional filtering using a wedge filter that has a smooth roll-off in direction.

[0052] As an example of bandpass filtering, a first order Butterworth bandpass filter or other suitable bandpass filter can be applied to each frame to remove spatial wavelengths that represent wave velocities outside a specified range. For example, the bandpass filter can be designed to remove spatial wavelengths that represent wave velocities outside a predetermined range of 0.5 to 1.0 m / s. Alternatively, the limits of the bandpass filter can be selected to follow a phase velocity curve, such as a phase velocity curve determined by a Fourier transform (e.g., 2DFT). For example, the limits can be selected to follow a phase velocity curve, such as a phase velocity curve determined by a Fourier transform (e.g., 2DFT). B (f)-c L ~c B (f)+c L So, the dispersion curve C B (f) ±C L can be selected by giving a range

[0053] This filtering process is performed in the frequency-wavenumber domain of the particle motion.

number

number

[0054] In the formula, B k (k Z ,k X ) is k Z -k X represents a two-dimensional filter of the region,

number

number

[0055] Spectral data for particular frequencies or frequency bands are then selected for processing, as shown in step 304. For example, the acquired spectral data

number

[0056] Next, the short-space two-dimensional ("2D") Fourier transform (F 2D ) is the windowed wavefield data, as shown in step 306.

number

number

[0057] In some implementations, the short-space 2DFT can be repeated for each pixel in the x and z directions. In some other implementations, the short-space 2DFT can also be repeated in a sparse manner every n-th or m-th pixel in the x and z directions, respectively. When the short-space 2DFT is performed in such a sparse manner, the missing pixels can be reconstructed using linear interpolation or other suitable steps. As an example, the linear interpolation step can be performed by applying a 2D convolution method.

[0058] In general, performing a short-space 2DFT on selected spectral data is

number

number

number

number

[0059] As an example, the window function may be a 2D cosine tapered window (ie, a Tukey window), which may be defined as follows:

number

[0060] During the ceremony,

number

number

[0061] where ξ and h denote the center and width of the window, respectively.

number

[0062] In general, increasing the spatial window size can decrease the average phase velocity. This behavior is caused by smoothing at the edges of the inclusions. When larger spatial windows are used, the phase velocity is smoothed significantly. Thus, a smaller spatial window size improves the spatial resolution of the wave velocity estimation. However, choosing a spatial window that is too small will result in a noisy image with a blurred form. Thus, the selection of the spatial window size can be chosen depending on the inclusion size to be evaluated. One approach is to gradually increase the spatial window size for image reconstruction. Using this approach, an optimal window size can be selected, which should be smaller than the inclusion under investigation to avoid excessive flattening and blurring of the edges.

[0063] As the spatial window slides along the spatial dimension, the windowed wave field domain

number

number

[0064] The resulting spectrum

number

[0065] In another implementation, a 1D Fourier transform can be applied along the spatial dimensions z and x separately, resulting in a set of 1D wavenumber spectra for a particular frequency f0, as follows:

number

[0066] It can be expressed in the z direction as:

number

[0067] In the case where, as in other examples of this disclosure, directional and / or bandpass filtering is applied to the wave data,

number

number

[0068] Next, the spatial distribution of the phase velocity of the wave versus frequency f is calculated, as shown in step 310. In an example where the wavenumber spectrum data is generated using a 2D Fourier transform,

number

number

[0069] During the ceremony,

number

number

[0070] and k (z:m) and k (x:m) is the argument calculated using:

number

[0071] Alternatively, in the example where the wavenumber spectrum data is generated using a discrete 1D Fourier transform, the spatial distribution of the wave phase velocity versus frequency f0 can be expressed as the wavenumber spectrum for all spatial locations as follows:

number

number

[0072] where the magnitude of the wavenumber is calculated as follows:

number

[0073] and in the formula:

number

number

[0074] Based on equations (21) or (24), a full field of view ("FOV") of a phase wave velocity image can be reconstructed, as shown at 312. As noted above, the frequency f0 can be expressed as a single frequency value, or as a frequency band f0 centered at f0. band , i.e., (f0-f b ,…,f0,…,f0+f b When using a frequency band, the phase velocity can be defined as

number

[0075] In the formula, f i (i=1,2,...,N) are frequency bands centered around f. Therefore, equation (28) can express the phase wave velocity image at each location over the selected frequency band.

[0076] These example methods have been described with respect to input of 3D wave data (e.g., two spatial dimensions, one time dimension). As noted above, the methods can also be extended for use with 4D wave data (e.g., three spatial dimensions, one time dimension). In these cases, the above process can be extended to include the y-axis (k y It can contain data along a single axis.

[0077] As a general matter, the phase velocity of viscoelastic materials (e.g., inclusions) varies with frequency and often plateaus at higher frequencies. For example, the phase velocity of inclusions may plateau above 850 Hz. Thus, for phase velocity imaging of inclusions or other viscoelastic materials, a wide bandwidth of frequencies may be advantageous. In most cases, a single acoustic radiation force push is used, which can provide waves of different frequencies depending on the medium. However, to obtain accurate measurements of wave velocity and other mechanical properties, it may be necessary to excite the waves at higher frequencies.

[0078] The phase velocity imaging techniques described in this disclosure can be used to reconstruct images at multiple frequencies to determine whether a wave velocity plateau has been reached over a defined frequency range. If this has not been reached, various techniques can be utilized to increase the frequency content of the propagating waves (e.g., shear, guided, or other mechanical waves). The acoustic radiation force push beam can be modified in its spatial distribution (e.g., by changing the number of elements used and / or focal position) or in its temporal application (e.g., by using repeated pushes at a particular frequency). To improve the accuracy of the measurements, feedback between image reconstruction and wave acquisition can be performed to optimize the bandwidth of the resulting wave. In this way, the phase velocity imaging techniques described in this disclosure can be used to assess the quality of the data and the reliability of the estimation of mechanical properties from a given data acquisition. Corrected acquisition settings can be generated in response to this analysis and used for subsequent acquisitions within a given imaging session.

[0079] Based on the local wave velocity image over a given frequency range, a mechanical property map can be generated. As an example, elasticity and viscosity maps can be generated from the local shear wave velocity image. For example, a KV rheological viscoelastic model can be used to generate the mechanical property map. The KV model consists of a dashpot and a spring arranged in parallel. The stress-strain relationship for the KV model can be expressed in the form:

number

[0080] In the formula, the stress σ is the product of the shear modulus μ1, the shear viscosity μ2, and the time derivative

number

number

[0081] In the formula, p represents density and ω is the angular frequency (ω=2πf).

[0082] Thus, by way of example, the elasticity and viscosity parameters can be estimated locally for each pixel by estimating the curve in equation (30) over the frequency domain, such as by using a nonlinear least-squares problem. For example, the following nonlinear least-squares problem can be solved:

number

[0083] Exemplary System 4, an example of a system 400 for generating a phase velocity image is shown, in accordance with some embodiments of the systems and methods described in this disclosure. As shown in FIG. 4, a computing device 450 can receive one or more types of data (e.g., ultrasound data, optical imaging data, magnetic resonance imaging data, mechanical wave data) from an image source 402, which can be an ultrasound imaging system, an optical imaging source, or an ultrasound image source, such as a magnetic resonance imaging system. In some embodiments, the computing device 450 can execute at least a portion of a phase velocity estimation system 404 to generate a phase velocity image from the data received from the image source 402.

[0084] Additionally or alternatively, in some embodiments, the computing device 450 can communicate information regarding the data received from the image source 402 to the server 452 via the communications network 454, and the server 452 can execute at least a portion of the phase velocity estimation system 404 to generate a phase velocity image from the data received from the image source 402. In such embodiments, the server 452 can return information indicative of an output of the phase velocity estimation system 404 to the computing device 450 (and / or any other suitable computing device) to generate a phase velocity image from the data received from the image source 402.

[0085] In some embodiments, the computing device 450 and / or the server 452 may be any suitable computing device or combination of devices, such as a desktop computer, a laptop computer, a smartphone, a tablet computer, a wearable computer, a server computer, a virtual machine run by a physical computing device, etc. The computing device 450 and / or the server 452 may also reconstruct an image from the data.

[0086] In some embodiments, the image source 402 can be any suitable source of image data (e.g., measurement data, images reconstructed from measurement data), such as an ultrasound imaging system, an optical imaging system, and an MRI system, another computing device (e.g., a server that stores image data), etc. In some embodiments, the image source 402 can be local to the computing device 450. For example, the image source 402 can be incorporated into the computing device 450 (e.g., the computing device 450 can be configured as part of a device for capturing, scanning, and / or storing images). As another example, the image source 402 can be connected to the computing device 450 by a cable, a direct wireless link, etc. Additionally or alternatively, in some embodiments, the image source 402 can be located locally and / or remotely from the computing device 450 and can communicate data to the computing device 450 (and / or the server 452) via a communications network (e.g., the communications network 454).

[0087] In some embodiments, the communication network 454 may be any suitable communication network or combination of communication networks. For example, the communication network 454 may include a Wi-Fi network (which may include one or more wireless routers, one or more switches, etc.), a peer-to-peer network (e.g., a Bluetooth network), a cellular network (e.g., a 3G network, a 4G network, etc., conforming to any suitable standard, such as CDMA, GSM, LTE, LTE Advanced, WiMAX, etc.), a wired network, etc. In some embodiments, the communication network 108 may be a local area network, a wide area network, a public network (e.g., the Internet), a private or semi-private network (e.g., a corporate or university intranet), any other suitable type of network, or any suitable combination of networks. Each of the communication links illustrated in FIG. 4 may be any suitable communication link or combination of communication links, such as a wired link, an optical fiber link, a Wi-Fi link, a Bluetooth link, a cellular link, etc.

[0088] 5, an example of hardware 500 that can be used to implement the image source 402, the computing device 450, and the server 454 according to some embodiments of the systems and methods described in this disclosure is shown. As shown in FIG. 5, in some embodiments, the computing device 450 can include a processor 502, a display 504, one or more inputs 506, one or more communication systems 508, and / or a memory 510. In some embodiments, the processor 502 can be any suitable hardware processor or combination of processors, such as a central processing unit ("CPU"), a graphics processing unit ("GPU"). In some embodiments, the display 504 can include any suitable display device, such as a computer monitor, a touch screen, a television, etc. In some embodiments, the input 506 can include any suitable input device and / or sensor that can be used to receive user input, such as a keyboard, a mouse, a touch screen, a microphone, etc.

[0089] In some embodiments, communications system 508 may include any suitable hardware, firmware, and / or software for communicating information over communications network 454 and / or any other suitable communications network. For example, communications system 508 may include one or more transceivers, one or more communications chips and / or chipsets, etc. In more particular examples, communications system 508 may include hardware, firmware, and / or software that may be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.

[0090] In some embodiments, memory 510 may include any suitable storage device(s) that may be used to store instructions, values, data, etc., and may be used, for example, by processor 502 to present content using display 504, communicate with server 452 via communication system(s) 508, etc. Memory 510 may include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 510 may include RAM, ROM, EEPROM, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 510 may have encoded thereon or otherwise stored therein a computer program for controlling the operation of computing device 450. In such embodiments, processor 502 may execute at least a portion of the computer program to present content (e.g., images, user interfaces, graphics, tables), receive content from server 452, transmit information to server 452, etc.

[0091] In some embodiments, the server 452 may include a processor 512, a display 514, one or more inputs 516, one or more communication systems 518, and / or memory 520. In some embodiments, the processor 512 may be any suitable hardware processor or combination of processors, such as a CPU, a GPU, etc. In some embodiments, the display 514 may include any suitable display device, such as a computer monitor, a touch screen, a television, etc. In some embodiments, the input 516 may include any suitable input device and / or sensor that may be used to receive user input, such as a keyboard, a mouse, a touch screen, a microphone, etc.

[0092] In some embodiments, communications system 518 may include any suitable hardware, firmware, and / or software for communicating information over communications network 454 and / or any other suitable communications network. For example, communications system 518 may include one or more transceivers, one or more communications chips, and / or chipsets, etc. In more particular examples, communications system 518 may include hardware, firmware, and / or software that may be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.

[0093] In some embodiments, memory 520 may include any public storage device(s) that may be used to store instructions, values, data, etc., and may be used by processor 512, for example, to present content using display 514, to communicate with one or more computing devices 450, etc. Memory 520 may include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 520 may include RAM, ROM, EEPROM, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 520 may have server programs encoded thereon for controlling operation of server 452. In such embodiments, processor 512 may execute at least a portion of the server programs to transmit information and / or content (e.g., data, images, user interfaces) to one or more computing devices 450, receive information and / or content from one or more computing devices 450, receive instructions from one or more devices (e.g., personal computers, laptop computers, tablet computers, smartphones), etc.

[0094] In some embodiments, the image source 402 may include a processor 522, one or more image acquisition systems 524, one or more communication systems 526, and / or memory 528. In some embodiments, the processor 522 may be any suitable hardware processor or combination of processors, such as a CPU, a GPU, etc. In some embodiments, the one or more image acquisition systems 524 are generally configured to acquire data, images, or both, and may include an ultrasound transducer, an optical imaging system, one or more radio frequency coils in an MRI system, etc. Additionally or alternatively, in some embodiments, the one or more image acquisition systems 524 may include any suitable hardware, firmware, and / or software for coupling to and / or controlling the operation of an ultrasound transducer, an optical imaging system, one or more radio frequency coils in an MRI system, etc. In some embodiments, one or more portions of the one or more image acquisition systems 524 may be removable and / or replaceable.

[0095] Although not shown, it should be noted that image source 402 may include any suitable inputs and / or outputs. For example, image source 402 may include input devices and / or sensors that may be used to receive user input, such as a keyboard, a mouse, a touch screen, a microphone, a track pad, a track ball, etc. As another example, image source 402 may include any suitable display device, such as a computer monitor, a touch screen, a television, etc., and one or more speakers.

[0096] In some embodiments, communications system 526 may include any suitable hardware, firmware, and / or software for communicating information to computing device 450 (and, in some embodiments, via communications network 454 and / or any other suitable communications network). For example, communications system 526 may include one or more transceivers, one or more communications chips and / or chipsets, etc. In more particular examples, communications system 526 may include hardware, firmware, and / or software that may be used to establish a wired connection using any suitable port and / or communications standard (e.g., VGA, DVI video, USB, RS-232, etc.), a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.

[0097] In some embodiments, memory 528 may include any suitable storage device(s) that may be used to store instructions, values, data, etc., which may be used by processor 522, for example, to control one or more image acquisition systems 524 and / or receive data from one or more image acquisition systems 524, to convert data into images, to present content (images, user interfaces, etc.) using a display, to communicate with one or more computing devices 450, etc. Memory 528 may include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 528 may include RAM, ROM, EEPROM, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 528 may have programs encoded thereon or otherwise stored therein for controlling the operation of image source 402. In such an embodiment, the processor 522 may execute at least a portion of a program to generate images, transmit information and / or content (e.g., data, images) to one or more computing devices 450, receive information and / or content from one or more computing devices 450, receive instructions from one or more devices (e.g., a personal computer, a laptop computer, a tablet computer, a smartphone, etc.), and the like.

[0098] 6 illustrates an example of an ultrasound system 600 capable of implementing the methods described in this disclosure. The ultrasound system 600 includes a transducer array 602 including a plurality of individually driven transducer elements 604. The transducer array 602 may include any suitable ultrasound transducer array, including a linear array, a curved array, a phased array, etc. Similarly, the transducer array 602 may include a 1D transducer, a 1.5D transducer, a 1.75D transducer, a 2D transducer, a 3D transducer, etc.

[0099] When energized by the transmitter 606, a given transducer element 604 generates a burst of ultrasonic energy. Ultrasonic energy reflected back to the transducer array 602 (e.g., echoes) from the object or subject under inspection is converted by each transducer element 604 into an electrical signal (e.g., echo signal) that can be individually applied to a receiver 608 via a set of switches 610. The transmitter 606, receiver 608, and switch 610 operate under the control of a controller 612, which may include one or more processors. As an example, the controller 612 may include a computer system.

[0100] The transmitter 606 can be programmed to transmit unfocused or focused ultrasound. In some configurations, the transmitter 606 can also be programmed to transmit diverging waves, spherical waves, cylindrical waves, plane waves, or combinations thereof. Additionally, the transmitter 606 can be programmed to transmit spatially or temporally encoded pulses.

[0101] The receiver 608 can be programmed to implement a detection sequence appropriate for the imaging task at hand. In some embodiments, the detection sequence can include one or more of line-by-line scanning, compound plane wave imaging, synthetic aperture imaging, and compound divergent beam imaging.

[0102] In some configurations, the transmitter 606 and receiver 608 can be programmed to implement a high frame rate, for example a frame rate associated with an acquisition pulse repetition frequency ("PRF") of at least 100 Hz. In some configurations, the ultrasound system 600 can sample and store at least 100 ensembles of echo signals in the time direction.

[0103] The controller 612 can be programmed to implement an imaging sequence as described in this disclosure or as otherwise known in the art. In some embodiments, the controller 612 receives user input that defines various factors used in implementing the imaging sequence.

[0104] A scan may be performed by setting the switches 610 to their transmit positions, thereby directing the transmitters 606 to momentarily turn on and energize the transducer elements 604 during a single transmit event in accordance with an imaging sequence. The switches 610 may then be set to their receive positions, and subsequent echo signals generated by the transducer elements 604 in response to one or more detected echoes are measured and applied to the receiver 608. The individual echo signals from the transducer elements 604 may be combined at the receiver 608 to generate a single echo signal.

[0105] The echo signals are communicated to a processing unit 614, which may be implemented by a hardware processor and memory, to process the echo signals or images generated from the echo signals. As an example, the processing unit 614 may generate a phase velocity image using methods described in this disclosure. The images generated from the echo signals by the processing unit 614 may be displayed on a display system 616.

[0106] In some embodiments, any suitable computer-readable medium may 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 transitory or non-transient. For example, non-transient computer-readable media may include media such as magnetic media (e.g., hard disks, floppy disks), optical media (e.g., compact disks, digital video disks, Blu-ray® disks), semiconductor media (e.g., random access memory ("RAM")), flash memory, electrically programmable read-only memory ("EPROM"), electrically erasable programmable read-only memory ("EEPROM")), any suitable media that is not transitory or lacks any semblance of persistence during transmission, and / or any suitable tangible media. As another example, a transitory computer-readable medium may include signals on a network, in wires, conductors, optical fibers, circuits, or any suitable media that is transitory and lacks any semblance of persistence during transmission, and / or any suitable intangible media.

[0107] The invention has been described in terms of one or more preferred embodiments, and it is to be understood that, except as expressly stated, many equivalents, alternatives, variations, and modifications are possible and are within the scope of the invention.

Claims

1. 1. A method for generating a phase velocity image from mechanical wave data acquired with an imaging system, the method comprising the steps of: (a) using a computer system to access mechanical wave data acquired by an imaging system; (b) estimating phase velocity data from the mechanical wave data using the computer system; and (c) generating a phase velocity image from the phase velocity data; Step (b) converting the mechanical wave data into frequency data by applying a Fourier transform to the mechanical wave data; selecting frequency components from the frequency data; windowing the frequency data associated with the frequency components using a sliding window to generate windowed wave field data; generating wavenumber spectrum data by Fourier transforming the windowed wavefield data; generating the phase velocity data as a spatial distribution of phase velocities estimated from the wavenumber spectrum data, thereby estimating the phase velocity data.

2. The method of claim 1 , wherein the mechanical wave data is converted to the frequency data by applying a one-dimensional Fourier transform along the time dimension to the mechanical wave data.

3. The mechanical wave data is applying a three-dimensional Fourier transform to the mechanical wave data to generate frequency-wavenumber data; filtering the frequency-wavenumber data to generate filtered frequency-wavenumber data; and applying a two-dimensional Fourier transform to the frequency-wavenumber data filtered along a first wavenumber dimension and a second wavenumber dimension to generate the frequency data.

4. The method of claim 3 , wherein filtering the frequency-wavenumber data comprises applying a directional filter to the frequency-wavenumber data to separate the frequency-wavenumber data into different directional components.

5. 4. The method of claim 3, wherein filtering the frequency-wavenumber data comprises applying a bandpass filter to the frequency-wavenumber data to remove spatial frequency data outside a predetermined spatial frequency range.

6. The method of claim 1 , wherein the frequency components consist of a single frequency.

7. The method of claim 1 , wherein the frequency components include a frequency band centered about a selected frequency.

8. The method of claim 1 , wherein a two-dimensional cosine tapered window is used to window the frequency data relating to the frequency components.

9. 2. The method of claim 1 , wherein the wavenumber spectral data is generated by applying a two-dimensional Fourier transform to the windowed wavefield data, the two-dimensional Fourier transform being applied along a first spatial dimension and a second spatial dimension.

10. 2. The method of claim 1 , wherein the wave number spectral data comprises first wave number spectral data and second wave number spectral data, the first wave number spectral data being generated by applying a first one-dimensional Fourier transform to the wave field data windowed along a first spatial dimension, and the second wave number spectral data being generated by applying a second one-dimensional Fourier transform to the wave field data windowed along a second spatial dimension.

11. The method of claim 1 , wherein the spatial distribution of the phase velocities is generated based in part on dividing the frequency components by the magnitude of the wavenumber spectrum data.

12. 2. The method of claim 1 , wherein step (b) comprises generating directionally filtered component data by directionally filtering the mechanical wave data into different directionally filtered components, and estimating phase velocity data for each directionally filtered component.

13. The method of claim 12 , wherein generating the phase velocity image comprises combining the phase velocity data estimated for each different directionally filtered component.

14. The method of claim 1 , wherein the mechanical wave data includes at least one of shear wave data or guided wave data.

15. The method of claim 1 , wherein the imaging system is one of an ultrasound system, a magnetic resonance imaging (MRI) system, or an optical imaging system.

16. The method of claim 1 , further comprising generating a mechanical property map from the phase velocity image.

17. The method of claim 16 , wherein the mechanical property map is generated based on fitting data in the phase velocity image to a model.

Images