Patent Description:
Elastography has been used in recent times to estimate biomechanical properties such as stiffness of a region of interest (ROI) in a patient, for example to assess the stiffness of a mass in a breast of abdomen as an aid in deciding if the mass is benign or cancerous or to determine other characteristics of the lesion. In principle, a force is applied to the ROI and the way the ROI deforms or moves in response, as measured by an imaging modality such as ultrasound or MRI, is used as an indication of stiffness. The known techniques include (a) quasistatic elastrography, in which the change in shape of a region of interest and its surroundings due to compression is measured, (b) acoustic radiation force impulse imaging (ARFI), in which the way a focused ultrasound beam pushes the ROI along the beam direction is used as a measure of stiffness, (c) shear wave elasticity imaging (SWEI), which is similar to ARFI but the ROI displacement normal to the beam direction is the measure of stiffness, and (d) supersonic shear imaging (SSI), in which acoustic radiation force pushing the ROI generates shear waves and the speed of those waves is a measure of stiffness. In magnetic resonance elastography, the speed of shear waves in the patient is used as a stiffness measure.

In the more detailed discussion of background below, and in the description of the new approaches, references are identified in parenthesis; full citations of the references are provided at the end of the specification. All cited referenced material are hereby incorporated by reference.

There has been a robust development of techniques to estimate and image the biomechanical properties of tissues (Parker et al. The applied stimulus can be quasi-static, transient, or continuous waves. Each one of these has unique mathematics and techniques for inverse solutions (Doyley, <NUM>), but all lie on a continuum of biomechanical responses (Parker et al. , <NUM>; Parker, <NUM>). Shear wave propagation has received significant attention, but an inherent challenge in many approaches is the presence of reflected waves from organ boundaries and internal inhomogeneities. These reflections are responsible for modal patterns (Parker and Lerner, <NUM>; Taylor et al. , <NUM>) in continuous wave applications and also for backwards travelling waves in transient experiments (Ringleb et al. Directional filtering can be used to eliminate some types of reflections (Pengfei et al. , <NUM>; McLaughlin and Renzi, <NUM>; Deffieux et al. , <NUM>; Manduca et al. , <NUM>; Engel and Bashford, <NUM>; Hah et al.

Continuous shear wave inversion approaches have been developed to estimate the unknown tissue stiffness. They include inversions of the Helmoltz equation in magnetic resonance elastography (MRE) (Van Houten et al. , <NUM>; Ringleb et al. , <NUM>; Romano et al. , <NUM>; Sinkus et al. , <NUM>; Oliphant et al. , <NUM>) and sonoelastography (Parker and Lerner, <NUM>; Fai et al. , <NUM>; Fu et al. Another class of estimations has been developed for underwater acoustics and geomechanics using random signals (Roux et al. , <NUM>) and has been extended to noise correlation measurements in soft tissues (Gallot et al. , <NUM>; Catheline et al. , <NUM>; Brum et al. These involve spatial coherence of noise functions measured at two points, and can be recast as Green's functions and time reversal solutions. An approach using a mechanical vibration source to produce multiple wave directions in a large organ such as the liver, has been developed by Tzschatzsch et al. (<NUM>; <NUM>; <NUM>). Using a probability approach, they characterize the shear wave speed by finding the minimum wavenumber as a function of direction. Elastographic systems and methods, using multiple shear waves, which interfere with each other are known from the documents <CIT>, <CIT>, <CIT>, <CIT>, <CIT> and <CIT>.

While the known techniques can be useful and beneficial, a need still remains for new approaches to estimating properties such as stiffness of a region of interest in a body that can provide improved performance and results. For example, there is a need for more effective ways to measure stiffness of regions of interest deeper in a patient, such as deep in the abdomen, and particularly in obese patients. There also is a need to reduce undesirable effects of acoustic energy interactions with tissue other than the region of interest. And, there is a need to generally improve performance and reduce the cost of equipment and its operation. This patent specification describes such new approaches.

The invention is defined in the independent claims <NUM> and <NUM>. This patent specification describes new approaches to measuring and estimating elastic properties of regions of interest or anatomical organ inside a body. The general framework involves using reverberant shear waves in bounded elastic media.

A non-limiting example of a system using a reverberant shear wave field in a body to measure viscoelastic properties of hidden regions of interest in the body comprises a plurality of vibration sources configured to produce shear waves at a selected frequency that interact with each other and with structures in the body to produce the desired reverberant shear wave field, an imaging device configured to measure motion of the region of interest in the body in a first selected direction in the presence of said reverberant shear wave field and to produce an estimate of the measured motion, an image processor configured to receive as inputs the selected frequency of the shear waves and said estimate of measured motion and to process the inputs with computer algorithms to provide an estimate of one or more viscoelastic properties of said region of interest in the body, a computer display configured to display said estimate of one or more viscoelastic properties of the region of interest, and a controller operatively coupled with said vibration sources, imaging device, and computer display to control their operation.

The imaging device can be an ultrasound scanner that includes an ultrasound transducer configured to provide a time sequence of ultrasound images of the region of interest and the estimate of measured motion. As another example, the imaging device can be a magnetic resonance imaging (MRI) machine configured to provide a time sequence of magnetic resonance images of the region of interest and the estimate of measured motion. Other examples of an imaging device include an x-ray imaging device that takes a rapid succession of x-ray images of the region of interest, and an Optical Coherence Tomography (OCT) device that similarly takes a rapid sequence of images. The vibration sources can be configured to produce shear waves that have substantially the same frequency, or are within the same range that in turn is within a wider range of, for example, <NUM>-<NUM> or more, such as <NUM>-<NUM>, and an even wider range such as <NUM>-<NUM>. In one example, <NUM>-<NUM> individual vibration sources can be used, but more than <NUM> can be used in other examples.

The estimate of measured motion includes a position in the selected direction of the region of interest at plural respective times and speed of the change in position in that direction. The image processor can be configured to calculate one or more estimates of the motion and speed as a function of auto-correlation of the positions of the ROI at respective times. The system can be further configured to measure motion and provide an estimate of the measured motion in one or two additional directions, in which case a map can be produced to show a spatial distribution of one or more viscoelastic parameters such as stiffness in a plane in 2D or a volume in 3D.

A detailed description of examples of preferred embodiments is provided below. While several embodiments are described, the new subject matter described in this patent specification is not limited to any one embodiment or combination of embodiments described herein, but instead encompasses numerous alternatives, modifications, and equivalents, as far as they fall within the scope of the claims. In addition, while numerous specific details are set forth in the following description to provide a thorough understanding, some embodiments can be practiced without some of these details. Moreover, for clarity and conciseness, certain technical material that is known in the related art has not been described in detail in order to avoid unnecessarily obscuring the new subject matter described herein. It should be clear that individual features of one or several of the specific embodiments described herein can be used in combination with features or other described embodiments. Further, like reference numbers and designations in the various drawings indicate like elements.

First, a discussion is presented of a theoretical basis for the new approaches to estimating stiffness of an internal ROI, and then examples of specific implementations are described.

<FIG> illustrates a geometry and a coordinate system that might help in appreciating the discussion below. The figure shows an imaging transducer <NUM> and an object <NUM> that has an isotropic random distribution of shear waves propagating through its interior, represented by the <MAT> direction vectors 104b of the individual plane waves that are distributed throughout the reverberant interior and in the region of interest 104a, consistent with a theory of reverberant shear wave fields. The x-axis of the coordinate system is aligned with the axial direction of the imaging transducer <NUM>, and it is assumed that the imaging system detects motion in the x-direction.

The complex pressure P̂ at a position ε in a reverberant chamber can be thought of as the superposition of plane waves incident from all directions (Pierce, <NUM>; Parker and Maye, <NUM>) <MAT> where the index q represents direction, nq are unit vectors uniformly distributed around <NUM>π solid angle, k and ω are the wave number and radial frequency of the plane waves, and P̂q are independent, identically distributed variables of random magnitude and phase. The corresponding velocity at a point is thereby given as <MAT> Where, from the plane-wave impedance relations, <MAT> where ρ is the media density, and c the speed of sound.

To calculate the autocorrelation function, the x component of velocity at some position ε within the tissue can be written as: <MAT> where êx is a unit vector in the x direction and <MAT> The summation on q is understood to be taken over <NUM>π solid angle.

Writing the correlation function definition then substituting equation (<NUM>), gives: <MAT> where E{ } represents an ensemble average and the asterisk represents conjugation. The product of the two series will include cross terms of the form: <MAT> But since the nxq and υ̂q are independent and the υ̂q are uncorrelated, this term vanishes. Thus: <MAT> Taking the real part of the equation (<NUM>), gives: <MAT> where, since the Vq are independent of the nqx and cosine terms, the mean squared value of the velocity is taken out from the curly braces. Since an ideal, diffuse field is assumed to be present in the reverberant chamber, then the ensemble or spatial averaging will assign equal weighting to all directions of incident sound. Thus, the average of the summation over discrete directions becomes the average over all directions of incident waves (Pierce, <NUM>; Cook et al. , <NUM>), around the polar coordinates of <FIG>: <MAT> Without loss of generality, vector Δε can be aligned with the z axis in <FIG>. Using spherical coordinates: <MAT> <MAT> and the differential solid angle is <MAT> so<MAT> Integrating first over φ and expanding the cosine term yields <MAT> where j<NUM>(x) is the spherical Bessel function of the first kind, of order <NUM>. This result is commensurate with the role of spherical Bessel functions in solutions to the Helmholtz equation via Fourier and Hankel transformations (Baddour, <NUM>). Also, equation (<NUM>) can be written in terms of trigonometric functions or Bessel functions of order <MAT> (Parker and Maye, <NUM>; Abramowitz and Stegun, <NUM>).

Now, switching to shear waves, the major difference is that the direction of propagation is perpendicular to the direction of displacement. Thus, if nq is taken as the direction of propagation, nqp is a perpendicular direction of shear displacement and velocity. Therefore, nq ·nqp = <NUM>.

To account for the perpendicular relation in the case of shear wave, <NUM>° or π/<NUM> is added to the angle formed by q and x (the detected direction). Thus, equation (<NUM>) becomes <MAT> and, following the same logical progression as before, equation (<NUM>) becomes <MAT>.

In the case where Δε is taken along the x -axis (the direction of the detected velocity), then nq · êxΔεx = Δεx sin θ cos φ , the argument in equation (<NUM>) becomes cos(ω<NUM>Δt - kΔεx sinθcosφ), and the integration results in <MAT>.

The two functions from equations (<NUM>) and (<NUM>) are shown in <FIG>, which illustrates that autocorrelation Bυυ (t = <NUM>, Δε) depends on the relative direction of Δε with respect to the detected component of the velocity vector. In this case, the detected direction of motion is taken as the x - axis. The thick line gives the autocorrelation with respect to x, the thin line shows the autocorrelation when Δε is taken in the perpendicular z - direction. The wavenumber k is set to unity for these curves.

The simplicity of equations (<NUM>) and (<NUM>), basically "sinc" and "jinc" spatial functions, is useful for practical implementations. An ultrasound or magnetic resonance imaging (MRI) scanner can track tissue motion within a region of interest. This generates a function v(x) along some region of interest (ROI). The tissue is subjected to multiple shear wave sources that are operating at a frequency typically in the range of <NUM> - <NUM>. The correlation function Bvv is calculated and fit to equation (<NUM>) or (<NUM>) to estimate the unknown parameter k. Local estimates of k are used to create a map, typically displayed in color, representing the shear wave speed and hence the stiffness of the tissue at different locations.

An efficient estimator for the unknown k in equation (<NUM>) is realized by examining the Fourier transform of the autocorrelation function: <MAT> where s is the spatial transform variable. This is a strictly bandlimited function with upper limit of spatial frequency set by k, the unknown wavenumber. The second moment m<NUM> of the transform is therefore similarly determined by k. From Bracewell, chapter <NUM>, page <NUM> (<NUM>): <MAT> Similarly, the Fourier transform of equation (<NUM>)'s spatial term is: <MAT> and the function is a real and even function of s. The second moment for this case is: <MAT> Furthermore, it is well known that the second moment of a transform is precisely related to the second derivative of the function at the origin (Bracewell, <NUM>). This can be approximated by a finite difference. Thus: <MAT> where k̂ is the estimate, C is a constant inversely dependent on Δx<NUM> , and the Δx lag and zero lag values of the real part of the autocorrelation at Δt = <NUM> from some segment of data are used. A similar expression applies to the estimate using Δz.

The new approaches described in this patent specification have been confirmed through numerical simulations discussed below.

Numerical simulations using finite elements analysis have been conducted using Abaqus/CAE version <NUM>-<NUM> (Dassault Systems, Vélizy-Villacoublay, France) in order to corroborate gelatin phantom experiments (shell-element analysis), and breast phantom (3D solid finite element analysis) described further below in this patent specification.

Figure <NUM> illustrates simulations where (a) shows a 3D mesh of a breast model, (b) shows a geometrical model of the breast shaped as a 2D profile of the mesh model in (a) with boundary conditions and loads indicated in Abaqus/CAE, (c) shows a finite shell-element mesh of (b) with material properties definition in the background and inclusion (diameter=<NUM>), and delimitation of the region on interest, and (d-f) shows a displacement field in the x-axis for the frequencies of operation: <NUM>, <NUM>, and <NUM>, respectively.

In Shell-element analysis, the profile of a 3D mesh model of a breast (Figure <NUM>(a)) is used to create a 2D geometrical solid model of a homogeneous breast with a hard inclusion simulating a tumor. The boundary conditions were set to have eight surface traction loads located in distinct parts of the breast model in order to produce shear displacement at the frequency of operation (Figure <NUM>(b)). The model was meshed using approximately <NUM>,<NUM>,<NUM> hybrid and quadratic shell elements (Figure <NUM>(c)). Power law frequency dependent viscoelastic material properties were chosen to represent tissue. The elastic modulus, density, and viscoelastic parameters assigned to the background and the inclusion are described in Table <NUM>. These parameters are in the range of gelatin phantom values extracted using the mechanical testing described below. The type of simulation was selected to be steady-state dynamic direct solution for a range of frequencies between <NUM> and <NUM>. This type of analysis calculates the 3D complex sinusoidal steady state solution of displacement and particle velocity of the breast when the defined loads introduce vibrations in the model.

Table <NUM>. Viscoelastic material parameters of the background and inclusion in the 2D finite element model. The power law frequency-dependent model is <MAT>, where g*(ω) is the Fourier transform of the non-dimensional shear relaxation function, <MAT> is a complex constant, a is a real constant, and f = ω/2π.

After the simulation was conducted, the complex values of particle velocity were stored for a posterior post-processing step. Figure <NUM>(d-f) shows reverberant vector fields within the background and inclusion for different frequencies.

In 3D solid finite element analysis, a 3D geometrical solid model of a homogeneous breast with a hard inclusion was created and meshed using approximately <NUM>,<NUM> hybrid and quadratic tetrahedral elements using the shape of the mesh model in Figure <NUM>(a). The material properties were chosen according to the specifications of the zerdine breast phantom (model <NUM>, CIRS Inc. , Norfolk, Virginia) assigned to the background and inclusion part as described in Table <NUM>.

Table <NUM>. Viscoelastic material parameters of the background and inclusion in the 3D finite element model. The power law frequency-dependent model is <MAT>, where g*(ω) is the Fourier transform of the non-dimensional shear relaxation function, <MAT> is a complex constant, a is a real constant, and f = ω/<NUM>π.

The type of simulation selected was a steady-state dynamic direct solution for two frequencies of operation <NUM> and <NUM>. The boundary conditions were set to be zero displacement in the sector that represents the chest wall. In addition, eight surface traction loads were located in different parts of the breast model in order to produce shear displacement at the frequency of operation. The complex values of particle velocity were stored for post-processing. A profile cut of the model shows reverberant vector fields within the background and inclusion for the frequencies of operation of <NUM> and <NUM> (<FIG>(a,b)). In this example, a 3D finite element solid model in Abaqus/CAE of a breast with magnitude displacement field of the complex sinusoidal steady-state solution within an interior plane was used. <FIG> illustrates a displacement field in the z-axis for the frequency of operation of <NUM>, and Figure <NUM>(b) is a similar illustration for <NUM>, in the background and inclusion (diameter=<NUM>) regions.

Experiments have been performed to verify certain aspects of the new approach described above.

<FIG> shows a setup using a zerdine breast phantom <NUM> in which a reverberant shear wave field is induced by two mechanical vibrators 502a, 502b and two miniature vibration sources 504a, 504b. An imaging ultrasound transducer <NUM> is a part of an ultrasound imaging system (not seen in the figure). A power amplifier (model <NUM>, Bruel & Kjaer, Naerum, Denmark), and a digital power amplifier (model LP-<NUM> A+, Lepai, Bukang, China) driven by a dual channel function generator (model AFG3022B, Tektronix, Beaverton, OR, USA) provided input signals to the two mechanical vibration sources (model <NUM>, Bruel & Kjaer, Naerum, Denmark), and the two miniature vibration sources (model NCM02-<NUM>-<NUM>-<NUM> JB, H2W, Linear Actuator, Santa Clara, CA, USA) vibrating at frequencies between <NUM> and <NUM> in contact with a gelatin-based and a zerdine breast phantom. A stiffer cylindrical inclusion <NUM> in diameter, embedded in an otherwise homogeneous background, was constructed following the procedure used by Hah et al. One phantom with gelatin (<NUM> Bloom Pork Gelatin, Gelatin Innovations Inc. , Schiller Park, IL, USA) concentrations of <NUM>% for the background and <NUM>% for the inclusion, was created by heating a mixture of gelatin, <NUM> of degassed water, <NUM> of Na-Cl and <NUM> of agar to <NUM>. The mixture was then cooled to approximately <NUM> and poured into a cubic mold (<NUM> × <NUM> × <NUM><NUM>) and was then allowed to rest at <NUM> overnight. Before the experiment, this phantom was taken out of its molds and left at room temperature for <NUM> hours. The size and shape of the zerdine phantom simulates a patient in the supine position, and it contains several solid masses that are at least two times stiffer than the background. Lesions range in size from <NUM> to <NUM> in diameter and are randomly positioned throughout the background. Additionally, a Verasonics ultrasound system (V-<NUM>, Verasonics Inc. , Kirkland, WA, USA), which enables high frame rate acquisition and a coherent plane wave compounding acquisition scheme, and a linear array ultrasound transducer <NUM> (Model L7-<NUM>, Verasonics Inc. , Redmond, WA, USA) were used to track the induced displacements using Loupas' estimator (Loupas et al. A <NUM>-D matrix of IQ data was stored for post-processing. A movie of average axial displacement is computed from the acquired <NUM>-D IQ data. In all the experiments the center frequency was <NUM> and the tracking pulse repetition frequency (PRF) is set to acquire at least <NUM> samples per cycle, i.e., PRF = <NUM> times the vibration frequency.

For the gelatin-based phantom materials, compression tests were performed on three cylindrical samples (approximately <NUM> in diameter and <NUM> in length) made with the same mixtures used to construct the gelatin-based media. A QT/<NUM> mechanical device (MTS Systems Co. , Eden Prairie, MN, USA) with a <NUM> N load cell was used to measure the stress-strain response. The compression rate was adjusted to <NUM>/s. These conventional mechanical measurements were considered the reference when assessing the elasticity properties of the cylindrical phantom.

Results from the numerical simulations and experiments are discussed below.

Shell-element analysis results: Complex-value displacement frames of the reverberation pattern within an anterior ROI containing the lesion were stored for analysis. <FIG> shows the real part of the displacement pattern for an excitation frequency of <NUM> in which the background and inclusion can be identified by the global separation size of the peaks and valleys of the interference. <FIG> shows a 2D shear wave speed (SWS) map obtained by applying the new approach to an autocorrelation window (<NUM> × <NUM><NUM>) which is moved to cover the entire ROI. Profiles in the x and y direction are taken from the 2D cross correlation, which are in good agreement with theory when compared to the equations (<NUM>) and (<NUM>) in <FIG>(c,d), respectively. Subsequently, the second moment approach described in equation (<NUM>) is applied to the correlation profiles taking into account the scale factor C obtained previously from simulations. The estimation of k̂ in both axes is averaged and used in the calculation of the shear wave speed by cs = ω/k̂ave.

<NUM>-D solid finite element analysis results: Complex-value displacement frames of the reverberation at <NUM> and <NUM> were obtained during simulations. <FIG> shows a profile of the real part of the displacement pattern for <NUM> excitation frequency. The effects of the viscosity (attenuation) and the boundary conditions influence the quality of the reverberation pattern. The same approach for the calculation of SWSM on the shell-element simulation data is applied for the 3D finite element case as shown in <FIG>(b,c). Average values of shear wave speed reports cs = <NUM> ± <NUM>/s in the background (ground truth cs = <NUM>/s) and cs = <NUM> ± <NUM>/s in the inclusion (ground truth cs = <NUM>/s). Some artifacts are seen near the outer boundaries where the assumptions of the isotropic reverberant fields may not hold.

Ultrasound experiments results when using the Gelatin-based phantom: <FIG> shows the reverberation displacements from the gelatin-based phantom at <NUM>. <FIG> presents the phase map obtained by the Fourier transform of the displacement pattern through time at each (axial-lateral) location. The B-mode and the SWS images (obtained by applying the new approach to a window <NUM> × <NUM><NUM> in size) are shown in <FIG>(c,d), respectively. A region of interest (<NUM> × <NUM><NUM>) was taken from the background and the inclusion in order to obtain a mean value of each region: <NUM> ± <NUM>/s, <NUM> ± <NUM>/s, respectively. The SWS result for both regions are in agreement with the elasticity properties obtained with mechanical measurements: <NUM> ± <NUM>/s, <NUM> ± <NUM>/s, respectively. Finally, <FIG> shows the SWS profile for a fixed depth corresponding to <FIG>.

Ultrasound experiments results when using the Zerdine breast phantom: <FIG> shows the reverberation results from the CIRS breast phantom. The displacement patterns using a vibration frequency of <NUM> is shown in <FIG> presents the phase map of the reverberation pattern at <NUM>; it can be noticed that there are two different wavelengths between <NUM> and <NUM> in the axial direction, which is the region where the inclusion is located (see B-mode image (<FIG>). The SWS maps, (obtained by applying our approach to a correlation window <NUM> × <NUM><NUM> in size) overlaying the B-mode image, are shown in <FIG>. A region of interest (<NUM> × <NUM><NUM>) was taken from the background and the inclusion in order to obtain a mean value for each region: <NUM> ± <NUM>/s, <NUM> ± <NUM>/s for <NUM>, respectively. The SWS result for the background is in agreement with the elasticity properties specified by the CIRS manufacturer: <NUM> ± <NUM>. For the SWS result at the inclusion, the phantom manufacturer only reported that the inclusion stiffness is, at least, two times higher than the background. Thus, the SWS result for the inclusion is in agreement with that information as well. Finally, <FIG> shows the SWS profile for a fixed depth corresponding to <FIG>.

<FIG> shows a SWS estimate summary, background and inclusion ROI for the gelatin-based and breast phantoms in portions (a) and (b) respectively. The width of the violin plots shows the probability density of the data at different values. The solid and dashed lines represent the mean and median of the data.

One example of a practical implementation of the new approach is described below, but it should be clear that this is one of several possible implementations and the new approach is not limited to this example.

As illustrated in block diagram form in <FIG>, a system <NUM> comprises a vibration source <NUM>, which in this example is in the form of individual vibration sources <NUM>-<NUM> through <NUM>-N, coupled with a body <NUM> and configured to produce shear waves that interact with each other and with structures in body <NUM> to produce a reverberant shear wave field in the body. In other examples, an integrated vibration source <NUM> can be used, such as an extended surface that moves differently at different points on the surface in a way that produces the required reverberant shear wave field in body <NUM>, and thereby functions similarly to an array of individual vibration sources. Arrows 104b designate randomly or nearly randomly distributed directions of shear waves in the reverberant shear wave field in body <NUM>. An imaging system <NUM> measures motion of a region of interest (ROI) 104a in body <NUM>, preferably in a single selected direction such as the beam direction of the imaging ultrasound transducer, in the presence of the reverberant shear wave field. An image processor <NUM> applies computer processing to the motion measurements from imaging transducer <NUM> and to parameters of the frequency of vibration sources <NUM> to produce estimates of one or more viscoelastic properties of ROI 104a such as stiffness properties. A computer display <NUM> displays the estimates and typically also includes an interface 110a through which a user communicates with the system. A controller <NUM> communicates with elements <NUM>, <NUM>, <NUM>, <NUM>, and 110a to monitor and control their operations.

An important aspect of system <NUM> is that it only needs to measure motion of ROI in a single direction, and that the shear waves in body <NUM> typically are in all or nearly all directions. Vibration source <NUM> can comprise several, typically <NUM> to <NUM> or more individual sources or points from which shear waves are emitted, arrayed around body <NUM> in a way that need not be precise so long as they can contribute to produce the specified reverberant shear wave field. The shear waves that source <NUM> produces need not be precisely directed. Individual vibration sources <NUM> can be used that operate independently of each other, or one or more integrated set or sets of sources <NUM> can be used. To generalize in the case of using individual vibration sources, vibration sources <NUM>-<NUM> through <NUM>-N are used, where N is a positive integer greater than <NUM> and preferably greater than <NUM>. In the case of an integrated vibration source, the source has <NUM> through N points or portions that produce respective shear waves that in turn produce the required reverberant field. Imaging system <NUM> can include an ultrasound transducer with a 1D or 2D array of transducer elements operating at a frequency suited to imaging the region of interest, such as <NUM> or another suitable frequency. As another example, imaging system <NUM> can use another imaging modality such as magnetic resonance, in which case the required imaging pulse sequence is simplified and is faster than for 2D or 3D MR imaging because only motion of ROI 104a in a single direction is required to be measured. As another example, imaging system <NUM> can employ Optical Coherence Tomography (OCT) or x-ray imaging. The images should be taken at the required time sequence frequency, for example two times or preferably more, such as five times or ten times, the highest frequency of interest of the reverberant shear wave field in the region of interest in the body. Image processor <NUM> can be a known ultrasound engine adapted to process the echoes from the ultrasound transducer into images that can be autocorrelated to derive displacement of ROI 104a and speed of the displacements in a selected direction, and thereafter the desired one or more viscoelastic properties of the ROI. In case imaging transducer <NUM> is a magnetic resonance scanner, image processor <NUM> can be the known computer facility of the MRI system, also adapted to provide measurements of displacement of ROI 104a and ultimately the desired one or more viscoelastic properties. In the case of OCT or x-ray imaging, image processor <NUM> can be the existing system's computer programmed to provide the images of the ROI from which can be derived the measurements of displacement and ultimately the desired viscoelastic properties. Computer display <NUM> can be the display that typically is a part of an ultrasound, MRI or OCT scanner or an x-ray imaging system. Controller <NUM> can be a separate device serving the indicated functions of can be implemented by suitably programming and interfacing an existing computer facility of an ultrasound or MRI or OCT scanner or x-ray imaging system.

In operation, a body <NUM> such as a patient's breast is placed on a platform <NUM> that is equipped with vibration a vibration source <NUM> such as individual sources <NUM>-<NUM> through <NUM>-N, for example <NUM>-<NUM> or more sources, that are coupled with breast <NUM> to create a reverberant shear wave field in it. The vibration sources can be on the platform, on a separate structure, or in a belt that can surround a patient's abdomen, or in some other integrated structure, so long as a sufficient area of body <NUM> is free to couple an imaging ultrasound transducer to the patient that can image the ROI and can induce the required reverberant filed body <NUM>. In the case MRI is used to image the ROI, suitable precautions need to be taken to account for the presence of a strong steady magnetic field and magnetic gradients in selecting and using the vibration sources and in their placement relative to the patient's body and MRI components such as RF coils. Other imaging modalities can be used instead of ultrasound and MRI, such as OCT (Optical Coherence Tomography) or x-ray imaging so long as they can produce the required time sequence of images of the ROI that show its motion in the presence of the reverberant shear wave field.

Main steps of the process are illustrated in <FIG>. In step <NUM>, a user such as a health professional selects a frequency of the shear waves, or a frequency band, and inputs the selection into controller <NUM> through interface 110a. The selected frequency preferably is in the range of <NUM>-<NUM>, and can be <NUM> as one example, but can be a different range, such as without limitation <NUM>-<NUM>, or <NUM>-<NUM>, for example when assessing bone in the patient may be of interest. Alternatively, the frequency or frequency range or ranges can be pre-set or otherwise communicated to the vibration source or sources so that the user need not input them for each patient study. Even if preset, the system can still allow the user to over-ride the preset and select a different vibration frequency or band or bands for a given patient study. Preferably the vibration source produces shear waves that have substantially the same frequency, or frequencies that are in a narrow band within a wider band such as <NUM>-<NUM> or another relatively wide range such as <NUM>-<NUM>, or <NUM>-<NUM>, but it is possible to produce shear waves that differ substantially in frequency and/or phase from each other, in which case the processing of the images of the ROI needs to take into account such substantial difference. In step <NUM>, the user images the ROI, for example by moving an ultrasound transducer that is a part of imaging system <NUM> over breast <NUM>, and observes the resulting images to select a position and an orientation of the transducer in which the ROI is within the image. For this purpose, the user can operate a commercially available ultrasound system such as those offered by GEHealthcare under the tradename LOGIQ, by Siemens (Acuson) offered under the tradename S3000, or another commercially available ultrasound scanner or a scanner purpose-designed specifically for system <NUM> of <FIG>. Such an ultrasound scanner includes an ultrasound transducer serving as an element of imaging system <NUM>, an ultrasound engine serving as element <NUM>, and a display and interface serving as elements <NUM> and 110a illustrated in <FIG>. In step <NUM>, the user turns on vibration source or sources <NUM> to produce the required reverberant shear wave field in breast <NUM> or at least at ROI 104a and surroundings. While the reverberant shear wave field is present in breast <NUM> or at least at ROI 104a and surroundings, in step <NUM> the ultrasound scanner continues to image ROI 104a to produce a rapid time sequence of images of the ROI. Image processor <NUM> receives from imaging system <NUM> the required images and in step <NUM> estimates at least one vector direction of the reverberant shear waves, Image processor <NUM>, which can be the ultrasound engine of the indicated commercially available ultrasound scanners, programmed as needed to carry out the operations described in this patent specification, receives as inputs the vibration frequency of sources <NUM> and the vector estimate in step <NUM>, and in step <NUM> computes the autocorrelation functions Bvxvx (Δt,Δεz) and Bvxvx (Δt,Δεx) identified in equations (<NUM>) and (<NUM>), respectively, using known autocorrelation algorithms to convert the measured displacements into the autocorrelation values. (Autocorrelations can be taken along one or more other Δε lines, instead of or in addition to taking autocorrelation along the Δεz and Δεx directions, for example along a diagonal line between the x and z axes or along some other direction that is not the direction of the imaging ultrasound beam or perpendicular to the beam. ) For this purpose, step <NUM> uses the displacements for increments Δt and Δεz from step <NUM>. In step <NUM>, image processor <NUM> computes an estimate of the wave number k and the shear wave speed c=w/k from the autocorrelation function per equation <NUM>. In step <NUM>, the system outputs to a computer display a map of the estimates shear wave speeds and/or viscoelastic parameters such as stiffness of the ROI or the organ of interest.

In principle, the calculation of the indicated autocorrelation functions Bvxvx (Δt,Δεz) and Bvxvx (Δt,Δεx) identified in equations (<NUM>) and (<NUM>) involves matching positions of the ROI or points of the ROI in successive ultrasound or MR or OCT images of the ROI and surroundings to determine how those positions have changed in shape or position from one point in time to another along the x-direction and/or along the z-direction in this example, and thus estimate motion in one or more directions of the ROI or points thereof. Bvxvx (Δt,Δεz) pertains to motion in the z-direction, for example perpendicular to the imaging beam from the ultrasound transducer, and Bvxvx (Δt,Δεx) pertains to motion in the x-direction, for example along the direction of the beam from the ultrasound transducer. Conventional processing, within the skill of an ordinary computer programmer, can be used to obtain the Δt, Δεz, and Δεx values from the succession of images of the ROI and use them to calculate the autocorrelation values and in turn use the calculated autocorrelation values in equation <NUM>. See Loupas, <NUM> for an example of autocorrelation processing. While autocorrelation along only one direction is in principle sufficient to derive the desired elastic property or properties of the ROI, such as stiffness or speed, in practical applications taking autocorrelation along two or even more directions may produce benefits such as reduction of noise effects resulting from redundancy of measurements. Once at least one of the autocorrelation values is available, equation (<NUM>) explains how to use it to derive the k (wave number) estimate. The coefficient C in equation (<NUM>) is empirically derived by tests with a known object such as a breast phantom where the parameters of equation (<NUM>) other than C are known or can be estimated, and solving for C. Stiffness is related to wave number k, as explained for example in Parker <NUM>, and can be derived for each point of interest in the ROI to estimate and display stiffness at a single point in the ROI, of as a 2D map of points in the ROI, or as a 3D map of points.

In review, there are several advantages to utilizing the framework of reverberant fields. First, the presence of reflections from boundaries and internal inhomogeneities is unavoidable in some common elastographic approaches, and these reflections plus the application of multiple sources and mode conversions can all contribute to the creation of the reverberant shear wave field. Once established, the characteristics of that field can be exploited to estimate the underlying shear wave phase velocity and/or viscoelastic properties. Secondly, the expected value of the autocorrelation function has been derived assuming only one vector component of detected velocity. This represents the simplest and most rapid data acquisition for both ultrasound and MRI, as additional transmit directions (in ultrasound) or additional phase encoding (in MRI) are required to determine additional vector components of shear wave velocity, and these are unnecessary in the framework developed herein. Thirdly, the need to verify a principal direction of wave propagation is eliminated in the reverberant framework as the underlying mathematics account for a superposition of waves. Fourthly, the need for explicit knowledge of boundary conditions or second derivatives that are essential in some inverse approaches (Doyley, <NUM>) are avoided in the reverberant approach.

Factors such as organ size, attenuation, frequency, and shear wave sources can cause degradation or deviation from the model. In addition, the performance of estimators of shear wave speed as a function of the same parameters can requires a more detailed assessment. Finally, multi-frequency versions of the new approach can be implemented to assess the frequency dependence of shear wave speed and hence the dispersion and viscoelastic properties.

Claim 1:
A system using a reverberant shear wave field in a body to estimate viscoelastic properties of hidden regions of interest in the body, comprising:
a vibration source (<NUM>) configured to produce plural shear waves at selected frequencies that interact with each other and with structures in the body to thereby produce a reverberant shear wave field in the body;
an imaging device (<NUM>) configured to measure motion of a region of interest in the body in a first selected direction in the presence of said reverberant shear wave field and to produce an estimate of the measured motion;
an image processor (<NUM>) configured to receive as inputs the selected frequencies of the shear waves and said estimate of measured motion, and to process the inputs with computer algorithms to provide an estimate of one or more viscoelastic properties of said region of interest in the body, using a framework of reverberant shear waves in bounded elastic media;
a computer display (<NUM><NUM>) configured to display said estimate of said one or more viscoelastic properties of the region of interest; and
a controller (<NUM>) operatively coupled with said vibration source, imaging device, image processor and computer display to control their operation.