Patent Description:
Non-invasive or non-destructive measurement of the mechanical properties of a medium is useful in a wide range of application. In particular, measuring the mechanical properties of tissues has important medical applications because it is related to tissue health state. For example, liver fibrosis is associated with increase of stiffness (shear modulus or shear elasticity) of liver tissue and thus measurement of liver stiffness can be used to non-invasively stage liver fibrosis. One way to non-invasively and non-destructively assess stiffness is using shear waves. As such, there has been increasing interest in creating and accurately measuring shear wave propagating in a medium.

Regardless of the particular system and resulting functionality being used or the underlying clinical information being sought, the use of shear waves in medical applications is increasing. As such, there is a need to provide more robust and efficient systems and methods for measuring or determining shear wave speed such in a manner appropriate for medical applications.

A non-invasive method for the quantitative determination of the mechanical properties of soft tissues in vivo known as Dynamic magnetic resonance elastography (MRE) is described in "Three-dimensional analysis of shear wave propagation observed by in vivo magnetic resonance elastography of the brain" by <NPL>.

In <CIT>, a system and method for measuring mechanical properties of a tissue using an ultrasound system is described. Ultrasound energy is applied to the tissue to produce shear waves that propagate in the tissue. Measurement data are then acquired by directing ultrasound detection pulses into the tissue. Information about the intensity field of the ultrasound energy used to produce the shear waves is obtained and used to produce a correction factor which is applied to the measurement data to correct the measurement data for errors arising from the geometry of the ultrasound energy used to produce the shear waves. From the corrected measurement data, mechanical properties of the tissue are calculated.

<CIT> describes methods for measuring mechanical properties of an object or subject under examination with an ultrasound system and using unfocused ultrasound energy are provided. Shear waves that propagate in the object or subject are produced by applying unfocused ultrasound energy to the object or subject, and measurement data is acquired by applying focused or unfocused ultrasound energy to at least one location in the object or subject at which shear waves are present Mechanical properties are then calculated from the acquired measurement data.

A method and system for imaging using virtual extended shear wave sources is described in <CIT>. Ultrasound energy is transmitted into tissue in a first direction to provide a virtual extended shear wave source. The virtual extended shear wave source generates an extended shear wave that propagates in a second direction substantially orthogonal to the first direction to cause movement in the first direction of tissue that is offset from the virtual extended shear wave source in the second direction.

A shear elasticity imaging technique, called comb-push ultrasound shear elastography (CUSE), in which multiple unfocused ultrasound beams arranged in a comb pattern (comb-push) are used to generate shear waves is described in <NPL>. A directional filter is then applied upon the shear wave field to extract the left-to-right (LR) and right-to-left (RL) propagating shear waves. A <NUM>-D shear wave speed map is then reconstructed by combining the LR and RL speed maps.

In "<NPL>, it is described that dynamic magnetic resonance elastography can visualize and measure propagating shear waves in tissue-like materials subjected to harmonic mechanical excitation. This allows the calculation of local values of material parameters such as shear modulus and attenuation. A spatio-temporal directional filter applied as a pre-processing step can separate interfering waves so they can be processed separately. Weighted combinations of inversions from such directionally separated data sets can significantly improve reconstructions of shear modulus and attenuation.

In accordance with one aspect of the present disclosure, a method of producing images of properties of an object is provided that includes producing a multi-directional shear wave field in the object. The method also includes using an imaging device, acquiring data about the multi-directional shear wave field in at least two spatial dimensions over at least one time instance and generating a set of component data by separating the data acquired into component data that contains isolated waves, each isolated wave propagating in a different one of a plurality of directions. The method further includes calculating from the component data for each wave contained in the component data at least a lateral shear wave speed component and an axial shear wave speed component, and producing a shear wave speed map for each propagation direction using the lateral and axial shear wave speed components. The method includes combining shear wave speed maps to produce at least one of a speed image and material property image for the object.

In accordance with another aspect of the present disclosure, a system is provided for producing images of properties of an object. The system includes an excitation system configured to produce a multi-directional shear wave field in the medium and a detection system configured to acquire data about the multi-directional shear wave field in at least two spatial dimensions over at least one time instance. The system also includes a processor configured to produce a multi-directional shear wave field in the object, and, using an imaging device, acquiring data about the multi-directional shear wave field in at least two spatial dimensions over at least one time instance. The processor is also configured to generate a set of component data by separating the data acquired into component data that contains isolated waves, each isolated wave propagating in a different one of a plurality of directions. The processor is further configured to calculate from the component data for each wave contained in the component data at least a lateral shear wave speed component and an axial shear wave speed component, and producing a shear wave speed map for each propagation direction using the lateral and axial shear wave speed components. The processor is also configured to combine shear wave speed maps to produce at least one of a speed image and material property image for the object.

The foregoing and other aspects and advantages of the invention will appear from the following description. In the description, reference is made to the accompanying drawings which form a part hereof, and in which there is shown by way of illustration a preferred embodiment of the invention. Such embodiment does not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims and herein for interpreting the scope of the invention.

Measuring the mechanical properties of tissues has important medical applications because it is related to tissue health state. For example, liver fibrosis is associated with increase of stiffness (shear modulus or shear elasticity) of liver tissue and thus measurement of liver stiffness can be used to non-invasively stage liver fibrosis. Propagation of shear waves is determined by a medium's mechanical properties. According to the Voigt model, shear wave speed cs in a medium relates to its shear modulus µ<NUM> and viscosity µ<NUM> by:
<MAT>.

Ultrasound can be used to generate shear waves remotely within the tissue for noninvasive elasticity imaging. Typically, a push ultrasound beam (focused or not focused) with long duration is used to produce a transient shear wave and pulse echo ultrasound is used to detect the propagation of the shear wave, as shown in <FIG>. Tissue particles move up and down due to the shear wave and this perturbation propagates outwards (see arrows in <FIG>) from the push beam at propagation speed cs. That is, in <FIG>, A, B, C, D, E are positions of shear wave detection by pulse-echo ultrasound. The vertical lines represent shear wave fronts, which are moving outwards from the push center. At the example shown here, the shear wave front (represented by the solid vertical line) is at position C. The dashed lines represent the shear wave fronts that have already passed (positions A and B) or to be arriving (positions D and E).

The time profile of shear wave motion at multiple positions detected by pulse-echo ultrasound along the shear wave propagation path can be used to calculate cs. For example, assuming the distance between position A and E is Δr and the delay between the arrival time of the shear wave at these two positions is Δt, then cs = Δr/Δt. The time delay Δt can be estimated by tracking the time instance of the shear wave peak at each position, or by finding the delay that gives the maximum cross correlation between the <NUM> shear wave time signals detected at each position. For the example shown in <FIG>, the shear wave fronts A-E are relatively uniform along depth direction z. Therefore, we only need to detect the shear wave along the x-direction to measure the shear wave propagation speed. In other words, a set of one-dimensional (1D) spatial data is required to correctly measure shear wave speed if the propagation direction and the detection direction are aligned.

Shear waves produced by ultrasound push beams are typically weak (micrometers), making shear wave detection, and therefore cs measurement, susceptible to noise (cardiac motion, breathing motion, body motion, ultrasound system noise, and the like). Therefore, shear wave measurements using ultrasound push beams have limited penetration. This can be problematic for applications requiring deeper penetration: for example, making shear wave measurements in livers of obese patients. Shear waves produced by mechanical vibration sources (external vibrator or internal cardiac motion) can be of much higher amplitude for more reliable measurements in deeper regions. For some applications, a multi-directional wave propagation field is desirable. To create a multi-directional vibration field one could use multiple small external vibrators that are activated in a continuous or transient manner. Additionally, a single large external vibrator could be used and driven with a continuous or transient excitation signal. Physiological motion such as that arising from the heart contraction, pressure waves in the vessels, or breathing can also be used as a source of wave motion. Each method has different advantages with respect to directionality of the waves, frequency characteristics, and motion amplitude. Shear waves thus produced typically come in multiple directions whose orientations are unknown. This can cause bias in shear wave speed measurements.

In the example shown in <FIG>, the shear wave is propagating at an oblique angle from wave front <NUM> to <NUM> to <NUM>. Assuming the wave front propagating from <NUM> to <NUM> (with a distance of a) within a time interval Δt, the actual shear wave speed cs is a/Δt. If shear wave speed is only measured along the x-direction as was the case in <FIG>, the apparent shear wave speed cs' will be b/Δt, which is higher than the real shear wave speed cs. Therefore, when the measurement direction is not aligned with shear wave propagation direction, the estimated shear wave speed will be biased high. This bias effect applies to the in-plane (within the transducer imaging plane) oblique shear waves shown in <FIG>, as well as out-of-plane (out of the transducer imaging plane) shear waves that travel through the two-dimensional (2D) ultrasound detection plane at an oblique angle.

Here two methods are provided to correct for this bias effect. The first method assumes the medium is homogenous. This can be used to study diffuse diseases where the mechanical properties are expected to change uniformly across the entire organ. Examples include liver fibrosis, lung fibrosis, and changes to the brain in patients with Alzheimer's disease.

Consider shear wave motion that can be measured in a 2D spatial plane in the x-z plane as shown in <FIG>. Referring to <FIG>, a process for analysis in a homogeneous media begins at process block <NUM> with acquiring a series of measurements repeatedly through time, t. Then, at process block <NUM>, the data is analyzed. For example, one way to analyze this spatiotemporal data, u(x,z,t), is to apply a Fourier transform along the two spatial dimensions and the temporal dimension to give U(kx,kz,f). This frequency domain representation will be referred to as k-f space.

For a given frequency, fc, the directionality of the wave motion can be examined by looking at the distribution of energy of |U(kx,kz,fc)|. An in-plane wave propagating at an arbitrary angle will show up as a peak in the k-f space with location (kx,kz,fc). To this end, at process block <NUM>, peaks are identified and used, at process block <NUM>, to determine a direction of wave movement. For example, a peak with a positive kX coordinate indicates a wave moving to right. A peak with a positive kZ coordinate indicates a wave moving downward. At frequency fc, a ray can be drawn from the origin of the kx-kz plane to the center of the peak, and the angle of that ray, θ = tan-<NUM>(kz/kx), indicates the wave propagation direction for the wave, and the radius of this component <MAT> is related to the shear wave speed cs by cs(fc) = fc/kr. At process block <NUM>, the location of peaks in k-f space are analyzed to determine similar waves, such as those with the same shear wave speed. For example, multiple waves, propagating at different oblique directions will show up as multiple peaks in the k-f space. At a given temporal frequency f, the peaks from multiple in-plane waves will have equal distance from the origin of the kx-kz plane and lie on a circle with kr = fc/cs because the medium is homogeneous and therefore shear wave speed should be identical for all propagation directions. Accordingly, at process block <NUM>, shear wave propagation in a homogenous medium can be reported without errors attributable to the above-described bias effect.

Referring to <FIG>, a case is illustrated where three waves, waves <NUM>, <NUM>, and <NUM>, are present. Wave <NUM> is on the circle <NUM> kr = fc/cs traveling upwards to the left. Wave <NUM> is also on the circle kr = fc/cs and has a larger energy and is traveling downward and slightly to the right. Wave <NUM> is an out-of-plane wave that has a high wave speed because it is located near the center of the kx-kz plane. Therefore, the speed of in-plane shear waves propagating at oblique angles can be correctly measured using 2D spatial data.

However, waves that are propagating obliquely to the plane of the measured motion (out-of-plane waves) are subject to a bias effect similar to that illustrated in <FIG>, and will be measured as propagating faster than the true shear wave speed of the medium. In k-f space, these out-of-plane waves are represented as peaks that lie within the circle <NUM> as kr = fc/cs. In a multi-directional wave field, there are multiple in-plane waves and multiple out-of-plane waves. At a given frequency f in the k-f space, this multi-directional wave field will show up as multiple peaks distributed on and within the circle kr = fc/cs. Therefore, the circle in the k-f space can be found with the largest radius (lowest wave speed) that sits on the outer edge of these multiple wave peaks, and use the radius of the circle to calculate the correct shear wave speed.

A method for finding the circle is provided below. The wave energy in all directions is integrated in a circle of radius kr in the k-f space using the relationships:
<MAT>
<MAT>
<MAT>.

The S(kr,fc) function has a sigmoid shape versus kr and has its steepest slope at the radial position that has the most energy added. This steepest slope point can be found by finding the maximum in dS(kr,fc)/dkr which is referred to as km. If there are sufficient in plane shear waves propagating at different directions, dS(kr,fc)/dkr would have a peak at the circle where all the in-plane shear waves sit (the lowest wave speed). This closely correlates with the shear wave speed of the medium as given by:
<MAT>.

Experiments were conducted in a homogenous elastic phantom by attaching multiple vibrators to the phantom walls. The vibrators provided an impulsive excitation and were activated in a random fashion to generate shear waves propagating in multiple directions. The process described above was used to evaluate the phase velocities at multiple frequencies, fc. Specifically, the k-f space distributions for the homogeneous phantom for three different shear wave frequencies,fc = <NUM>, <NUM>,<NUM> was studied. To this end, <FIG> shows the results from using Eq. (<NUM>) on the data acquired. The circles on the graphs depict the location of the maximum slope and the value of km used for calculating the phase velocity with Eq. (<NUM>). The integrated magnitude continues to increase with kr, even when kr is larger than km. This is because the background in <FIG> has non-zero values. The phase velocity dispersion for different frequencies is shown in <FIG> and the viscoelastic material properties were estimated as µ<NUM> = <NUM> kPa and µ<NUM> = <NUM> Pa·s by fitting the data to the Voigt model in Eq. (<NUM>). These values of µ<NUM> and µ<NUM> are close to results obtained by independent validation measurements in this phantom. Although 2D spatial data is used as an example above, the method can be extended to 3D spatial data.

Referring to <FIG>, for an inhomogeneous medium, a different approach can be used to obtain the spatial distribution of shear wave speed or shear modulus. If the wave motion field is multi-directional, that is, waves are traveling in many different directions, the waves propagating in different directions can be isolated by acquiring data at process block <NUM>. Optionally, at process block <NUM>, a curl operation can be used to remove compressional wave while keeping shear waves. At process block <NUM>, directional filters are applied. The directional filters may be applied in the k-f space by multiplying the directional filter response with the data. An example may include eight directional filters separated by <NUM> degrees. Reference is made to the following, "<NPL>.

To reduce the effects of the out-of-plane shear waves that are measured as waves with a high biased speed, we can additionally filter these out by filtering out the low spatial frequency (k) values for the propagating waves. This can be incorporated into the directional filters by setting a lower limit on the values of kl for each frequency fc such that speeds above c = fc/kl (i.e., kr < kl) are eliminated. Similarly, an upper limit of ku can be set for each frequency fc such that waves with speed below c = fc/ku (i.e., kr > ku) are eliminated. This lower wave speed limit can be used to remove false "wave motions" caused by body motion or other unwanted interference during shear wave data acquisitions.

For example, the shear wave speed of human liver (from normal to cirrhosis) should be in the range of <NUM>-<NUM>/s. In this case, the lower and upper limit of shear wave speed can be set to <NUM> and <NUM>/s to suppress interfering waves with propagation speed out of this range so that the final shear wave speed estimation is more reliable. For a diffuse disease such as liver fibrosis, the above-disclosed "k-f space method for homogeneous media" can be used to obtain an initial estimate of the shear wave speed of the medium, which can then be used to set the lower and upper shear wave speed limits. For example, in a particular patient, the shear wave speed estimated by "k-f space method for homogeneous media" is <NUM>/s. Then the upper speed limit can be set for this particular patient to <NUM>/s to improve the rejection of out-of-plane waves and produce a 2D image with less bias. A 2D image will allow the calculation of variation of shear wave speed estimation within the liver as an indication of measurement reliability, and therefore is still valuable even when imaging a homogenous medium. Instead of a lower and higher wave speed limit, fixed thresholds kl and ku can be used for all frequencies fc such that kr < kl or kr > ku are eliminated.

At very low frequency fc, the resolution of kr may not be sufficiently small to allow proper elimination of waves based on wave speed limits. One solution is to set a lower temporal frequency limit fl and eliminate all waves with frequency lower than fi. A smooth ramping profile can be used instead of a unit step function to reduce ringing effects during this process. By way of example, if the shear waves produced by the external vibrations are mainly above <NUM>, then eliminating all waves with frequency lower than <NUM> will remove unwanted low frequency motions while keeping the useful shear waves intact. Similarly, an upper temporal frequency limit fu can be used to remove high frequency noise for more stable results. The speed or frequency limits imposed by hard threshold may create jump discontinuities associated with undesirable Gibbs ringing artifacts. Therefore, a "soft" threshold with smooth transition can be used instead for these limits.

Referring again to <FIG>, after applying the directional filters at process block <NUM>, an inverse Fourier transform may be applied at process block <NUM> to bring the data back to the spatiotemporal domain for analysis to, at process block <NUM>, obtain shear wave speed or shear modulus. At process block <NUM>, the properties of the shear wave, including shear wave speed or shear modulus, are reported. As will be described below, this report may include shear wave speed maps and/or information from combine shear wave speed maps from multiple directions.

For example, the shear wave speed can be estimated using a time-of-flight (TOF) method (for example,<NPL>. ) or using normalized cross-correlation of the wave motion (for example, <NPL>. McLaughlin and <NPL>. These are one-dimensional methods, but in the imaging plane we can measure the shear wave speed in two-dimensions (2D), as shown in <FIG>.

In particular, a normalized cross-correlation can be applied to both x and z directions so that a lateral shear wave speed VX and an axial shear wave speed VZ can be obtained. By way of example, let the shear wave signal detected at pixel a, b, and c be Sa(t), Sb(t), and Sc(t), where t is time. Let the time delay estimated by cross-correlation be tab between Sa(t) and Sb(t), and tac between Sa(t) and Sc(t). Let distance between pixel a and c be Lac, and distance between pixel a and b be Lab. Then VX = Lac/ tac, and Vz = Lab/ tab. In the triangle denoted by apex a, b, and c, the true shear wave speed V can be calculated by the formula:
<MAT> or
<MAT>.

Equation (<NUM>) is more stable than Eq. (<NUM>) when either tac or tab is zero (if wave propagation direction is aligned with axis x or z). Note that the 2D vector shear wave speed calculation given by Eqs. (<NUM>) and (<NUM>) does not require a priori knowledge of the direction of shear wave propagation, which is difficult to know in practice.

Two methods were developed to increase the robustness of 2D vector shear wave speed calculation while preserving the spatial resolution. First, an algorithm used in numerical differentiation calculation was adapted to local shear wave speed calculation. Conventional local shear wave speed measurement techniques as introduced in <NPL>. , cross-correlate two shear waveforms (waveform of particle displacement or velocity versus time) from two imaging pixels that are a fixed distance apart, as shown in <FIG>, to estimate the shear wave speed of the center pixel. A more robust approach, as shown in <FIG>, cross-correlates multiple pairs of shear waveforms that are a shorter distance apart and produces multiple local shear wave speed estimates. The final shear wave speed at the center pixel is given by weighted summing these estimates by their correlation coefficients. This algorithm can be implemented along both x and z directions to obtain VX and VZ, as shown in <FIG>. With respect to <FIG>, all pixels within the gray window as in <FIG> may be used to get estimates of Vx and Vz. The triangles indicate the spatial locations of estimated Vx. The rectangles indicate the spatial locations of estimated Vz. The gradient shading indicates the distance weighting: higher weights are assigned to estimates that are closer to the center pixel (indicated by darker gray).

A total number of D shear wave speed maps are produced from D directional filters with different directions, using the 2D vector calculation method described above. The final shear wave speed map can be obtained from combining the D speed maps. In particular, <FIG> shows <NUM> speed maps (D = <NUM>) obtained in a test phantom with a hard inclusion in the middle, after applying the <NUM> directional filters. Multiple miniature vibrators were attached to the surface of the phantom to produce the complex wave field used in this experiment. <FIG> shows the final shear wave speed map obtained by weighted sum of the data in <FIG>. The weighting used here is the normalized cross-correlation coefficient (NCC) of each pixel in each directional image in <FIG>. Shear wave energy at each image pixel (the sum of the square of the motion signal over time) in each directional image is also a good quality control factor, and therefore can be used in combination with the NCC to calculate the weighting for each pixel. One can set up thresholds for shear wave energy and NCC so that shear wave speed estimates from pixels with low shear wave energy and poor cross-correlation will not contribute to the final map. Another option of combining multiple directional shear wave images is to take the minimum or median value of all images to produce the final shear wave map. This will suppress the bias caused by out of plane waves and compressional waves, and will work especially well for a homogenous medium.

After directional filtering, Eqs. (<NUM>) and (<NUM>) calculate the shear wave speed from two orthogonal directions. Shear wave data in real applications have noise and thus errors in time delay estimation in either x or z direction will enter the final shear wave estimation through Eqs. (<NUM>) and (<NUM>). In the presence of noise, it is desirable to measure shear wave propagation in multiple directions to improve robustness of final shear wave speed estimation. Referring to <FIG>, the shear wave front propagates from position <NUM> to <NUM> to <NUM> with a distance interval r. Assuming the time it takes for the shear wave to propagate from <NUM> to <NUM> is τ. Propagation from position <NUM> to <NUM> will also take time τ. The time delay between the shear wave signal at a (center of circle) and b (on the circle with radius r) is τab = τ · sin (θ), where θ is the angle determined by the location of b as shown in <FIG>. In reality, the measured delay is τab(θ) = τ · sin(θ) + n(θ), where n (θ) is noise. n (θ) at different angles should be independent to each other. Therefore, multiple delays τ (θ) at multiple angles can be measured and fit to a model τ · sin (θ) to obtain a more robust estimate of τ in the presence of noise n (θ). In practical situations, the propagation direction of the shear wave in the x-z coordinate is typically unknown. Therefore, the model used for fitting should be:
<MAT>
where φ is a constant phase offset. Both τ and φ can be determined during the data fitting process. And the shear wave speed is calculated by:
<MAT>.

When ultrasound detection of shear wave signal is performed on a Cartesian spatial grid, the delay measured between shear waves detected at <NUM> pixels needs to be scaled by the distance between these two pixels. Referring to <FIG>, the black dots denote the pixel location where shear waves are detected. Pixel a locates at the origin, and pixel b locates at (x, z). The distance R between a and b is therefore R = <MAT>, and the angle θ = arctangent (z/x). Assuming the delay measured between pixel a and b is Dab, the delay normalized by R is therefore τab = Dab/R. Using this approach, τab (θ) can be measured along multiple angles when other unlabeled black dots in <FIG> are selected as pixel b. The measured τab (θ) is then fit with the model in Eq. (<NUM>) to estimate τ, which can be used to calculated the final shear wave velocity with Eq. (<NUM>) by setting r = <NUM>, because the delays measured here have already been normalized by distance. <FIG> shows an example of data obtained in a test phantom with nominal shear wave speed of <NUM>/s. Dots are time delays (normalized by distance) measured at different angles. The solid line is the fitting of the data to a sine model, which gives an amplitude of <NUM>/m corresponding to a shear wave speed of <NUM>/s. The above model fitting approach uses the delay in multiple directions to calculate the final shear wave speed. Similar to the relationship between Eqs. (<NUM>) and (<NUM>), appearant wave speeds instead of delays in multiple directions can also be used to calculate the final shear wave speed. The approach shown in <FIG> can be combined with the model fitting approach to make the method more robust. In addition, each pixel will have a R<NUM> value indicating the quality of the model fitting with delay through Eq. (<NUM>) or with speed. Therefore, the R<NUM> value can be used to control the weighting of the shear wave speed maps for different directions to produce a final shear wave speed map.

Other than time-of-flight method and normalized cross-correlation method, phase lag method ("<NPL>) can also be used to estimate frequency dependent wave propagation speed (phase speed) using spatiotemporal data after directional filtering. In the phase lag method, a Fourier transform, Kalman filtering, or other appropriate methods are performed on the time signal at each pixel to calculate the phase of the shear wave at multiple frequencies. The shear wave speed at a given frequency f is then estimated from phase lag at frequency f of at least two pixels along the shear wave propagation direction. The phase lag method can resolve wave propagation speed at multiple frequencies (dispersion), and can be extended to calculate phase speed for waves propagating from unknown directions using Eq. (<NUM>). When time delays in Eqs. (<NUM>) through (<NUM>) are used instead, the time delay t can be computed from the phase lag p and the frequencyfof the shear wave: <MAT>. Although 2D spatial data is used as an example above, the method can be extended to 3D spatial data.

The above disclosed method needs all pixels of the spatiotemporal data to have the same time grid. This requirement can be met when "flash imaging" is used for pulse echo detection of shear waves. For traditional ultrasound scanners where the 2D data are acquired in a line-by-line or zone-by-zone manner, pixels at different ultrasound A-lines are sampled at different time grid. In such situation, the time signal at each pixel can be interpolated to a higher sampling rate such that the post-interpolation time samples are aligned on the same time grid for different pixels. One example of a system and method for interpolation and alignment is described in co-pending <CIT> and having publication number<CIT>.

Systems and methods of efficiently producing shear waves from external vibration are also disclosed here. For example, referring to <FIG>, a piece of furniture, such as a bed <NUM> or a chair <NUM> or other structure having a rigid surface or surfaces <NUM>, <NUM> can be used to vibrate the subject or medium (not shown) engaged with the rigid surface <NUM>, <NUM>. The rigid surface <NUM>, <NUM> can be flat or curved to optimize vibration transmission into the body. A vibration source <NUM>, <NUM> is coupled to the rigid surface <NUM>, <NUM>. As will be described, the vibration source <NUM>, <NUM> may include a motor with off-center weight or other appropriate sources of vibrations, such as pneumatic or acoustically delivered vibrations, that delivers vibrations to a subject or medium through the rigid surface <NUM>, <NUM>. In this regard, the vibration source may include a connection <NUM>, <NUM> that connects to a driving or power source. To this end, as will be described, depending upon the design of the vibration source <NUM>, <NUM> the connection <NUM>, <NUM> may be an electrical connection or may be a pneumatic or acoustic connection or other means of receiving driving energy to operate the vibration source <NUM>, <NUM>.

Cushions <NUM> can be added between the rigid surface <NUM>, <NUM> and the body to improve patient comfort. In addition to traditional cushion materials such as foams and rubbers, liquid or air sealed inside a flexible membrane can also serve as a deformable cushion that can conform to the body surface to maximize contact area for more efficient vibration delivery. Because the air or liquid is sealed in a closed space, vibration from the rigid surface can transmit through the air or liquid cushion efficiently into the body. The air or liquid cushion can be integrated with the rigid surface <NUM>, <NUM>.

It is contemplated that the vibration source <NUM>, <NUM> may include a large vibration source in some designs, such as a motor carrying an offset load. In other configurations, instead of one large source of vibrations, multiple vibration sources may be positioned at different locations along the rigid surface <NUM>, <NUM> to produce shear waves in the body. Alternatively, vibration sources can be secured directly to a body surface of the subject or medium, such as by elastic strap, adhesive membrane, or other appropriate means. Further still, the vibration sources can be embedded in cushions <NUM> for the patient to lie on, or embedded in a vest or other clothing for the patient to wear.

As mentioned DC motors with off-centered weight is one example of a system that may be used. In some cases, the DC motors may be as small as those used cell phone vibrators. Alternatively, larger electromagnetic actuators can be used to produce stronger shear waves in deeper regions from the body surface. Also, a pneumatic or other driver used by Magnetic Resonance Elastography (MRE) can be used for this purpose. Examples of such MRE driver systems may be found in <CIT>; <CIT>; <CIT>; <CIT>; and <CIT> and Application Nos. <CIT> and <CIT>. When ultrasound is used for shear wave detection, electromagnetic interference is not of concern. Thus, electromagnetic devices can be used as vibration sources.

Regardless of the specific means of generating or driving the vibrations, there are some basic configurations that can be used to select generation or driving sources for the vibration sources. For example, as shown in <FIG>, one type of actuator <NUM> has a "push" motion when activated. Such actuators <NUM> can be pressed against the body surface, for example, using elastic straps or other appropriate means. Upon activation, an actuation body <NUM> pushes, by way of a shaft <NUM>, on a contact plate <NUM> to produce a small "punch" to the body surface to generate shear waves inside the body. In another configuration, referring to <FIG>, an actuator <NUM> may utilize a "pull" motion when activated. In this configuration, the actuator is inside the actuation body <NUM> is connected through a shaft <NUM> to the contact plate <NUM>. However, a spring or other appropriate biasing mechanism <NUM> may be inserted between the actuator body <NUM> and the contact plate <NUM> to maintain separation of these two parts. The contact plate <NUM> is fixed to a case <NUM> of the actuation body <NUM> and the actuator is free to move and slide along the shaft <NUM>. Upon activation, the actuator body <NUM> accelerates and moves towards the contact plate <NUM>. The inertial mass of the actuator body <NUM> will produce an impact to the patient body for shear wave generation when it strikes the contact plate <NUM>. When the actuator is not activated, the bias mechanism <NUM> will separate the actuator body <NUM> and the contact plate <NUM> to provide an adequate distance for acceleration of the actuator body <NUM> for the next impact.

As discussed, multiple sources of vibration may be used. When multiple sources of vibrations are used to produce shear waves in a patient study, electric current drawn from a common source driving these vibrators can be very high, requiring a high power source. To solve this issue, a control circuit may be used to connect the driving source and the multiple vibrators. The control circuit connects the source with each of the multiple vibrators only for a brief period of time. The "ON" time of each vibrator can be purposely misaligned such that at any given time instant, the source is only powering one or a few vibrators. By way of example, assume one source is driving <NUM> vibrators. Each vibrator is turned on for <NUM> milliseconds (ms) to produce shear waves in the patient body. The control circuit can sequentially turn on and off each vibrator such that vibrator <NUM> is turned ON during <NUM>-<NUM>, vibrator <NUM> is turned ON during <NUM>-<NUM>, vibrator <NUM> is turned ON during <NUM>-<NUM>,. , vibrator <NUM> is turned ON during <NUM>-<NUM>. In this manner, the source only needs to drive one vibrator at a time, thus reducing the current and power requirement of the source.

An air or underwater loudspeaker can also be used to introduce vibration in body. The frequency of sound emitted by the loudspeaker can be adjusted to the resonant frequency of the body part where vibrations are to be introduced. For example, resonance may help introducing vibration in liver for shear wave elasticity measurements. In this case, the loudspeaker can be embedded within an examination bed with the active surface facing up and roughly level with the bed surface. The patient can lie on the bed facing upwards with the upper back positioned on top of the loud speaker. The opening of the bed for embedding the loud speaker should be large enough to allow sufficient transmission of acoustic energy from the loud speaker into the body. The opening should also not be too large so that it is completely covered and preferably sealed by the patient's back when the patient lies on the bed. Frequency of the sound emitted by the load speaker can be adjusted to match the resonant frequency of the chest cavity of the patient, for example at <NUM>. Alternatively, a chirp signal (for example, from <NUM> to <NUM>) can be emitted into the body for shear wave production. Then vibration from the rib cage and lung will propagate with the body to produce multi-directional shear waves in liver.

The methods disclosed here assume a plane wave propagating in a homogenous medium. In practice, this assumption can be considered valid locally. For example, within a <NUM> by <NUM> area, the shear wave can be considered as plane wave and the medium can be considered homogeneous.

Although shear waves is used as an example in this teaching, the methods disclosed here can be used to measure speed of other waves such as compressional waves. And shear waves can be detected by methods other than ultrasound, such as magnetic resonance imaging (MRI) or optical systems. For all detection methods, wave motion component (usually a 3D vector) in the x direction, y direction, z direction, or their combination when available, can be used with the above method for shear wave speed measurements.

Referring now to <FIG>, an example of an ultrasound imaging system <NUM> that may be used with the present invention is illustrated. It will be appreciated, however, that other suitable ultrasound systems can also be used to implement the present invention. The ultrasound imaging system <NUM> includes a transducer array <NUM> that includes a plurality of separately driven transducer elements <NUM>. When energized by a transmitter <NUM>, each transducer element <NUM> produces a burst of ultrasonic energy. The ultrasonic energy reflected back to the transducer array <NUM> from the object or subject under study is converted to an electrical signal by each transducer element <NUM> and applied separately to a receiver <NUM> through a set of switches <NUM>. The transmitter <NUM>, receiver <NUM>, and switches <NUM> are operated under the control of a digital controller <NUM> responsive to the commands input by a human operator. A complete scan is performed by acquiring a series of echo signals in which the switches <NUM> are set to their transmit position, thereby directing the transmitter <NUM> to be turned on momentarily to energize each transducer element <NUM>. The switches <NUM> are then set to their receive position and the subsequent echo signals produced by each transducer element <NUM> are measured and applied to the receiver <NUM>. The separate echo signals from each transducer element <NUM> are combined in the receiver <NUM> to produce a single echo signal that is employed to produce a line in an image, for example, on a display system <NUM>.

The transmitter <NUM> drives the transducer array <NUM> such that an ultrasonic beam is produced, and which is directed substantially perpendicular to the front surface of the transducer array <NUM>. To focus this ultrasonic beam at a range, R , from the transducer array <NUM>, a subgroup of the transducer elements <NUM> are energized to produce the ultrasonic beam and the pulsing of the inner transducer elements <NUM> in this subgroup are delayed relative to the outer transducer elements <NUM>, as shown at <NUM>. An ultrasonic beam focused at a point, P , results from the interference of the separate wavelets produced by the subgroup of transducer elements <NUM>. The time delays determine the depth of focus, or range, R, which is typically changed during a scan when a two-dimensional image is to be performed. The same time delay pattern is used when receiving the echo signals, resulting in dynamic focusing of the echo signals received by the subgroup of transducer elements <NUM>. In this manner, a single scan line in the image is formed.

To generate the next scan line, the subgroup of transducer elements <NUM> to be energized are shifted one transducer element <NUM> position along the length of the transducer array <NUM> and another scan line is acquired. As indicated at <NUM>, the focal point, P , of the ultrasonic beam is thereby shifted along the length of the transducer <NUM> by repeatedly shifting the location of the energized subgroup of transducer elements <NUM>.

Referring particularly to <FIG>, the transmitter <NUM> includes a set of channel pulse code memories, which are indicated collectively at <NUM>. In general, the number of pulse code memories <NUM> is equal to the number of transducer elements <NUM> in the transducer <NUM>. These pulse code memories are also referred to as transmission channels for this reason. Each pulse code memory <NUM> is typically a <NUM> × N bit memory that stores a bit pattern <NUM> that determines the frequency of the ultrasonic pulse <NUM> that is to be produced. This bit pattern <NUM> may be read out of each pulse code memory <NUM> by a master clock and applied to a driver <NUM> that amplifies the signal to a power level suitable for driving the transducer <NUM>. In the example shown in <FIG>, the bit pattern is a sequence of four "<NUM>" bits alternated with four "<NUM>" bits to produce a five megahertz ultrasonic pulse <NUM>. The transducer elements <NUM> to which these ultrasonic pulses <NUM> are applied respond by producing ultrasonic energy. If all of the available bits are used, a pulse with a narrow bandwidth centered on the carrier frequency will be emitted.

Referring particularly to <FIG>, the receiver <NUM> is comprised of three sections: a time-gain control section <NUM>, a beam forming section <NUM>, and a mid-processor section <NUM>. The time-gain control section <NUM> includes an amplifier <NUM> for each receiver channel in the receiver <NUM>, and a time-gain control circuit <NUM>. The input of each amplifier <NUM> is connected to a respective one of the transducer elements <NUM> to receive and amplify the echo signal that is receives from the respective transducer element <NUM>. The amount of amplification provided by the amplifiers <NUM> is controlled through a control line <NUM> that is driven by the time-gain control circuit <NUM>. As the depth, or range, R, of the echo signal increases, its amplitude is diminished. As a result, unless the echo signal emanating from more distant reflectors is amplified more than the echo signal from nearby reflectors, the brightness of the image diminishes rapidly as a function of range, R. This amplification is controlled by a user who manually sets time-gain control potentiometers <NUM> to values that provide a relatively uniform brightness over the entire range of the sector scan. The time interval over which the echo signal is acquired determines the range from which it emanates, and this time interval is divided into, for example, eight segments by the time-gain control circuit <NUM>. The settings of the time-gain control potentiometers <NUM> are employed to set the gain of the amplifiers <NUM> during each of the respective time intervals so that the received echo signal is amplified in ever increasing amounts over the acquisition time interval.

The beam forming section <NUM> of the receiver <NUM> includes a plurality of separate receiver channels <NUM>. As will be explained in more detail below, each receiver channel <NUM> receives an analog echo signal from one of the amplifiers <NUM> at an input <NUM>, and produces a stream of digitized output values on an in-phase, I, bus <NUM> and a quadrature, Q, bus <NUM>. Each of these I and Q values represents a sample of the echo signal envelope at a specific range, R. These samples have been delayed in the manner described above such that when they are summed with the I and Q samples from each of the other receiver channels <NUM> at summing points <NUM> and <NUM>, they indicate the magnitude and phase of the echo signal reflected from a point, P , located at range, R, on the steered beam, θ.

The mid-processor section <NUM> receives beam samples from the summing points <NUM> and <NUM>. The I and Q values of each beam sample may be, for example, a <NUM>-bit digital number that represents the in-phase, I, and quadrature, Q, components of the magnitude of the echo signal from a point (R,θ). The mid-processor <NUM> can perform a variety of calculations on these beam samples, the choice of which is determined by the type of imaging application at task.

Claim 1:
A method of producing images of properties of an object, comprising steps of:
a) producing a multi-directional shear wave field in the object;
b) using an imaging device, acquiring (<NUM>) data about the multi-directional shear wave field in at least two spatial dimensions over at least one time instance;
c) generating a set of component data by separating (<NUM>) the data acquired in step b) into component data that contains isolated waves, each isolated wave propagating in a different one of a plurality of directions;
d) calculating (<NUM>) from the component data for each wave contained in the component data, at least a lateral shear wave speed component and an axial shear wave speed component;
e) producing a shear wave speed map for each propagation direction using the lateral and axial shear wave speed components; and
f) combining shear wave speed maps to produce at least one of a speed image and material property image for the object.