EURASIP Journal on Applied Signal Processing 2005:8, 1212–1220 c ○ 2005 Hindawi Publishing Corporation Beamforming Scheme for 2D Displacement Estimation in Ultrasound Imaging

We propose a beamforming scheme for ultrasound imaging leading to the generation of two sets of images, one with oscillations only in the axial direction and one with oscillations only in the lateral direction. Applied to tissue elasticity imaging, this leads to the development of a specific displacement estimation technique that is capable of accurate estimation of two components of the displacement. The mean standard deviation for the axial displacement estimates is 0.0219 times the wavelength of the axial oscillations, and for the lateral estimates, it is equal to 0.0164 times the wavelength of the lateral oscillations. The method is presented and its feasibility is clearly established by a simulation work.


INTRODUCTION
In medical ultrasound imaging, beamforming can have different aims. In systems generating 3D volumes, sparse synthetic aperture beamforming can be used to maintain a frame rate comparable with existing 2D scanners [1]. Other beamforming techniques can be aimed at estimation of new parameters, like those associated with tissue elasticity imaging [2]. More generally beamforming enables to control many aspects of the image formation like the depth of field or the diffraction of the transmitted beam, and so forth [3,4]. In this paper, beamforming is introduced in the field of tissue elasticity imaging with ultrasound.
Estimation of elasticity of biological soft tissue with ultrasound deals with the mapping of any parameter characterising the elastic properties of the medium. Examples are Poisson's ratio [5] or Young's modulus. The latter has been shown to be highly affected by various pathological conditions [6] and consequently, elasticity imaging carries potential for diagnosing these diseases.
The basic principle involved in tissue elasticity imaging with ultrasound or "elastography" is to acquire ultrasound images of a medium in two different states. The first image is used as a reference. Then the investigated medium is compressed, and a second image is acquired. As the medium typically exhibits variations in stiffness, not all of it will have the same deformation due to the external compression. It is this difference of deformation that the different elastography methods try to recover in order to highlight the elastic properties of the medium.
The first published approaches attempted to estimate the axial displacement or strain field in the medium [7,8]. Unfortunately, knowledge of this component of the strain field is not enough to recover Young's modulus in the medium. Thus, additional components of the displacement or the strain field have to be estimated. Classical 2D block matching or speckle tracking methods [9] have been studied by others. In order to track more precisely the acoustic signature, some approaches have calculated interpolated RF signals located between adjacent received RF signals [5]. However these methods all feature low lateral resolution.
The estimation of the displacement along the axial direction has been shown to be very accurate due to the oscillations that are naturally present in the RF signals. This paper proposes a method that provides-in the lateral direction of the image-the same oscillating characteristics as in the axial direction. The method uses a newly introduced beamforming scheme [10]. This scheme can provide two different point spread functions (PSF). A classical PSF, only oscillating in the axial direction, and a new PSF only oscillating in the lateral direction. The latter is obtained by heterodyning demodulation [11] of a PSF oscillating in both directions, a so-called double oscillating PSF (DOPSF). The beamformer producing DOPSF has been introduced in blood flow estimation [10,12]. The use of a PSF oscillating only in the lateral direction produces "lateral RF signals" showing the same oscillating characteristics as in the axial RF signals allowing high-resolution 2D parameter estimation. The introduction of two beamformers in the field of tissue elasticity imaging is a new approach. It leads to the development of a specific displacement estimation for the axial and lateral directions of the image. This paper shows the feasibility of accurate estimation of both axial and lateral displacement fields by successive estimation and compensation steps.
The paper will first focus on the beamforming method that yields the expected PSF. Next the image formation model is developed followed by a presentation of the axial and lateral estimation. Finally numerical simulation results are presented followed by a conclusion on the work.

BEAMFORMING METHOD
In our approach two PSFs are used. In this part the beamformers used to create each of them are described. The problem is considered in the image plane only.
In linear acoustic systems, the point spread function can be separated [13] such that one term corresponds to the emission process and the other one to the receive process: Ultrasound probe with N a active elements where h e (x, z) can be interpreted as the spatial distribution of the emitted acoustic energy, and h r (x, z) represents the spatial distribution of the collected energy. ⊗ z denotes the convolution over the spatial axial variable z. Lateral oscillations can be introduced in the PSF in different ways. The most convenient is obtained by beamforming the signals received on the aperture by all active elements during the receive process.
With this approach the emitted beam should have the smallest possible effect on the shape of the PSF in order to control the global PSF only during the receive process [10]. It is quite complicated to give an exact analytical solution for the beamformer that can produce a specified pulse-echo PSF in a broadband case. An interesting approach to beamformer design is given in [3,4]. Nevertheless, it is possible to design the beamformer using a simpler approach without degrading the results too much [10]. Therefore, in the present paper, the beamformer is designed using continuous wave (CW) emission and Fraunhoffer approximation. Fraunhoffer approximation is obtained by using quadratic focusing [14].
The elements of the ultrasound probe will be modelled as point sources. These sources emit a spherical sinusoidal pressure wave in the medium. The pressure at depth z and discrete lateral position x j due to the distribution of sources on the aperture can be expressed as where g(x i , x j , z) is a Green's function in infinite space [4] which corresponds to the propagation function from the ith pressure source at discrete lateral position x i on the aperture to the point of interest located at depth z and discrete lateral position x j , as illustrated in Figure 1. w(x i ) is an apodization applied to the ith source. N a is the number of active elements on the aperture. The complex expression of the emitted pressure in front of the transducer can be written as where r i j = z 2 + (x i − x j ) 2 is the distance between the ith source and the point of interest where the pressure P(x j , z) is calculated. λ is the wavelength of the emitted wave. In order to reach the condition of Fraunhoffer approximation at the depth z, the difference in path d xi with is introduced in (3), where z 2 + x 2 i is the path from the ith element to the point of interest and z is the path from the center of operture to the point of interest. This corresponds to quadratic focusing and leads to Next, Fresnel approximation [14] is used to replace r i j and d xi in (5), by their first-order binomial expansion given by leading to the following approximation of (5): Denoting the phase terms outside the sum with C(x j , z) gives Aside coefficient C(x j , z), the second term in (8) can be interpreted as a discrete Fourier transform of the apodization function w(x i ). This result is equivalent to Fraunhoffer approximation [14]. From this observation it can be deduced that the apodization function leading to the expected pressure profile, is the inverse discrete Fourier transform (IDFT) of this profile with a division of the frequency variable x j by λz: These expressions were deduced for the emitted beam, but according to the reciprocity theorem, (9) can also be applied to obtain the receive apodization function. Equation (8) shows that the frequency variable involved in the Fourier transform relation is scaled with respect to λ and z. This means that the apodization function has to be adapted dynamically during the receive process. This leads to a receive apodization function depending on the depth, w(x i , z), in order to maintain a sensitivity profile which is constant with depth.
The receive beamforming scheme has to produce two oscillating PSFs, a classical one, oscillating only in the axial direction, and another one, oscillating only in the lateral direction.
The two PSFs are obtained from the same array of received signals for each active element during the receive process. This array is denoted by b k (x i , z). It is a function of lateral position on the array x i and depth z, the latter corresponding to time. Indice k indicates the kth ultrasound emission.
The PSF oscillating only in the lateral direction is obtained by axial demodulation of a DOPSF. The apodization yielding this DOPSF was obtained by inverse Fourier transform of the profile of the expected PSF. The RF image of the DOPSF is obtained by juxtaposition of the beamformed sig- where w 1 (x i , z) is the apodization matrix, d xi is the delay matrix, x k is the kth RF line in the image of the DOPSF corresponding to the kth ultrasound emission, h(x, z), and z varies from zero to the maximum depth of investigation. The DOPSF can be described in space as a separable cosine multiplication, where f x and f z denote the lateral and axial spacial frequency of the oscillations, respectively: It is possible to suppress the axial oscillations by combining this PSF, with H z {h(x, z)} and H x {h(x, z)} which are the Hilbert transforms in the axial and lateral directions, respectively, as explained in [11] and recalled here: The "over line" indicates complex conjugation. Equation (14) shows a complex PSF. Its real part gives the expected PSF without axial oscillations whereas its imaginary part is its lateral Hilbert transform. It is also important to notice that the spatial lateral frequency of the PSF oscillating only in the lateral direction is twice that of the DOPSF. A second point spread function is obtained with a classical beamformer. The same quadratic dynamic focusing will be applied. But this time, an apodization window, w 2 (x i ), constant with depth, will be used: Here, h z (x, z) is the usual point spread function of an ultrasound scanner, with only axial oscillations.

IMAGE FORMATION MODEL
In this part a general model of the image formation is first presented. Using this model, it is shown that an ultrasound image of a displaced medium can be considered as a displaced version of the image acquired before the medium was displaced. Then the specific images for the displacement estimation obtained with the two beamformers are presented.

General model
The ultrasound RF image, composed by the juxtaposition of received beamformed RF signals is the spatial convolution between the PSF and the scatterer distribution that composes the medium to be imaged. Thus the RF image before displacement (reference image) is where h(x, z) and d(x, z) represent the PSF and the distribution of scatterer strengths, respectively. Assuming a uniform scatterer strength this can be written as where N s is the total number of scatterers and (x i , z i ) are the coordinates of the ith scatterer.
Let the global displacement vector applied to the scatterers be ∆ = (∆ x , ∆ z ). The distribution of scatterers after displacement becomes The image after displacement (strain image) is given by which means that It is clear in this expression that the strain image is a displaced version of the reference image. As a consequence the displacement in the tissue can be estimated by estimating the displacement that occurs between two successive images.

Specific images
Four images have to be derived. Two reference images which are obtained by beamforming the same set of raw signals received from the medium in relaxed state. And two images which are obtained by beamforming the second set of raw signals, received from the medium after it has been displaced. The two previously defined PSFs are going to be used to achieve the two images of each state of the medium. Thus is the reference image obtained with the PSF oscillating only in the axial direction (z), is the strain image obtained with the PSF oscillating only in the axial direction (z), is the reference image obtained with the PSF oscillating only in the lateral direction (x), and is the strain image obtained with the PSF oscillating only in the lateral direction (x). R{h x (x, z)} denotes the real part of h x (x, z).

DISPLACEMENT ESTIMATION
A displacement estimation algorithm is developed based on local time delay estimation between ultrasound RF signals [8]. Images are simulated using (21)-(24). An estimation of the local displacement is calculated over the entire image. This is done through a sliding window. After each estimation, the reference window is shifted. The corresponding window on the delayed signal is also shifted by a distance taking into account the estimated displacement.
This estimation can be done for the axial and for the lateral directions separately. But there can be some problems of decorrelation due to a displacement in both directions of the image. To overcome this problem, an estimation loop has been setup with estimation and compensation steps. The displacement estimated for one direction is compensated before estimating in the other direction. The steps of estimation and compensation can be repeated as an iterative process until a certain criterion of quality of the estimation has been met.

Axial estimation and compensation
The axial displacement is obtained with the first set of images obtained with the PSF oscillating only in the axial direction. The columns of the images defined in (21) and (22) are the signals considered. An axial cross-correlation function is defined: between the reference signal and the Hilbert transform of the strain signal. The axial autocorrelation of the reference signal is also defined: where M is the number of samples in one window of an axial signal. The local axial displacement is estimated between the corresponding signals in two consecutive images [8]. A ma-trix∆ z (x, z) is estimated iteratively: where ω z is the axial angular frequency and the superscript i in∆ i z (x, z) indicates the iteration number. Before estimating the lateral displacement, the axial displacement that has been estimated is compensated in the strain images. To do this, each pixel in the strain image is displaced with respect to the displacement estimated for its position. This leads to an axially compensated image. The compensation is made for both sets of images, the one obtained with the PSF oscillating only in the axial direction and the one obtained with the PSF oscillating only in the lateral direction:

Lateral estimation and compensation
Now the two images, r x (x, z) and s xac (x, z), obtained with the PSF oscillating only in the lateral direction are considered. The lateral signals, composed by the juxtaposition of samples from the same depth, are considered. Again, consider the lateral cross-correlation function between the reference signal and the Hilbert transform of the strain signal: and an autocorrelation function of the reference signal: where N is the number of samples in one window of a lateral signal.
The lateral displacement matrix∆ x (x, z) is defined iteratively [8] aŝ where ω x is the lateral angular frequency, and the superscript i, in∆ i x (x, z), indicates the iteration number. The lateral compensation is made by shifting laterally each pixel with respect to the displacement estimated for its position. This results in an axially and laterally compensated image If the estimated displacements are correct, the doublecompensated images are identical to the initial ones. If not, the loop can be iterated one more time. The images considered during the second iteration are the reference images and the axially and laterally compensated images (32).

NUMERICAL SIMULATION RESULTS
The method described has been tested with computer simulations. The ultrasound images have been calculated with the Field II program developed by Jensen [15,16]. A real ultrasound probe, a B-K Medical type 8560 probe for use with a B&K Medical 3535 ultrasound scanner, has been simulated. Parameters are given in Table 1.  Figure 2) for h x The desired lateral profile of the DOPSF has to contain some oscillations and has to be limited in space. As a consequence, it has been decided that the pressure profile for the DOPSF should correspond to the following analytical expression: where rect(x j /L) is a rectangular window defined as The apodization matrix is calculated using (9), that is, the inverse discrete Fourier transform of the desired pressure profile given in (33). This apodization function is shown in Figure 2. The inverse Fourier transform of the product of a sinusoïd and a rectangular window is a sinc function convolved with two delta functions. This is in accordance with the profile of our apodization matrix. The position of the delta functions on the active part of the probe, which corresponds to the position of the maximum peak of the sinc functions, is directly related to the frequency of the lateral oscillations [10]: where x 0 is equal to half of the distance between the two peaks.
For the example in Figure 2, at a depth of 15 mm, the distance separating the two peaks is equal to 16 elements, which is equal to 6.34 mm. This leads to x 0 equal to 3.17 mm and 1/ f x equal to 1 mm for the DOPSF lateral wavelength, with an axial wavelength of 0.22 mm. From the DOPSF, it is possible to create a PSF oscillating in the lateral direction only by taking the real part of expression (14). This PSF is seen in Figure 3. Its lateral frequency should be equal to twice that of the DOPSF. The simulated PSF is in accordance with the expected lateral wavelength of 0.5 mm.
The typical axially only oscillating point spread function is obtained as expressed in (15) and can be seen in Figure 4.
To validate our method, a homogeneous numeric phantom with known displacement was simulated. It was composed of a spatial uniform random distribution of scatterers located in a plane 35 mm wide by 10 mm deep. The Poisson ratio was 0.49, which is a typical value for biological tissues. The problem was considered only in the image plane. The phantom was located 10 -20 mm from the transducer. It was compressed axially by 2%. This lead to a lateral strain of 0.98%. The exact displacement distribution and the estimated one can be seen in Figure 5. Notice how the displacement is more and more important the deeper the depth.
With the displacement scheme chosen, the axial displacement component is supposed to be constant for a particular depth. To characterize the axial estimation, the mean value of the estimate and its standard deviation have been calculated for different depths. This is reported in Figure 6.
For the lateral estimation, the same representation can be done. Indeed for a particular lateral position, the lateral displacement is constant with respect to the depth. The mean value of the estimation and its standard deviation are calculated and reported in Figure 7.   Figures 6 and 7 show that the estimated profile is extremely close to the true one. The mean standard deviation for the axial estimate is equal to 0.0219 λ z and for the lateral estimate it is equal to 0.0164 λ x . The estimation shows that the axial precision is in the same order of magnitude as the lateral one, when expressed in numbers of axial or lateral wavelengths. Because the axial wavelength is shorter than the lateral wavelength, the result is still better for the axial estimate than for the lateral estimate in absolute terms. To improve on this, the frequency of the lateral oscillations should be increased. Theoretically there is no limit for this, except the width of the active part of the probe. This is explained by (35). If the probe was larger, a higher frequency could be reached for the lateral oscillations, and the results of the estimation in mm could be more similar for both directions.

CONCLUSION
In this paper, a new approach to displacement estimation by ultrasound for tissue elasticity imaging has been presented. This approach is based on a beamformer designed with the CW theory. The raw signals coming from all active elements on the aperture are processed to create lateral oscillations in the RF ultrasound image. This leads to the development of   a specific estimation loop for axial and lateral displacement estimation. This loop is based on successive displacement estimation/compensation steps on two sets of images. The feasibility of the method has clearly been demonstrated with results obtained with simulation studies. Indeed low bias and standard deviation are obtained for both directions of estimation. Moreover the lack of precision in the lateral direction has been highly improved with next-to-similar performances in the lateral and axial directions.