WO2007110669A1 - Image data processing systems - Google Patents

Abstract

This invention generally relates to methods, apparatus and computer program code for processing data captured from imaging systems, in particular pulse-echo imaging systems such as ultrasonic imaging systems, in order to determine deformation data for an imaged object. A method of processing at least one-dimensional image data captured by an imaging technique to determine deformation data defining an at least one-dimensional deformation in an imaged object, said deformation data comprising at least one data pair, said data pair comprising displacement data and corresponding displacement location data, the method comprising: inputting first and second sets of said image data, a said set of image data comprising imaging signal data including at least one of signal magnitude data and signal phase data for an imaged region of said object; positioning a window at a position on said first set of image data and determining a corresponding position for said window in said second set of image data; determining displacement data defining a displacement estimate from said positions of said window in said first and second image data; and determining an estimated location of said displacement estimate responsive to said imaging signal data within said window for at least one of said first and second sets of image data; whereby said at least one data pair comprises said displacement data and said estimated location of said displacement estimate. Figure 5b

Description

Image Data Processing Systems

This invention generally relates to methods, apparatus and computer program code for processing data captured from imaging systems, in particular pulse-echo imaging systems such as ultrasonic imaging systems, in order to determine deformation data for an imaged object.

Ultrasonic strain imaging is usually based on displacement estimates computed using finite-length sections of the RF ultrasound signal. Amplitude variations in the ultrasound cause a perturbation in the location at which the displacement estimate is valid, and if this goes uncorrected, it is an important source of estimation noise, which is amplified when the displacement field is converted into a strain image. We will describe a correction method, tested on phase and correlation coefficient strain imaging. Results indicate that the new correction yields a substantial reduction in estimation noise. Background prior art can be found in US 6,520,913, which shows time delay (strain) calculated on a uniform grid (Figure 1), US 2005/165309, and US 2003/0200036.

Ultrasonic elasticity imaging spans a broad range of techniques that process ultrasound signals to extract information relating to tissue's mechanical properties. A majority of these techniques require high quality displacement tracking at the first stage of signal processing. Examples include quasistatic compression imaging, axial shear wave imaging and acoustic radiation force imaging in both quasistatic/impulsive and dynamic forms. The principal alternative, sonoelasticity imaging, employs Doppler velocity estimation in mechanically vibrated tissues. This is a practical technique, although the images it yields are relatively difficult to interpret. Displacement-based imaging systems have been investigated for a wide range of diagnostic purposes, spanning screening for soft tissue tumours, monitoring of atherosclerosis, assessment of skin pathologies, and examination of cardiac disease among other applications. The simplest form of meaningful visualisation is the strain image. This may be extended by analysing strain image sequences to infer material property estimates such as elastic and viscoelastic moduli.

The cornerstone of elasticity imaging is displacement tracking. Consider a pair of ultrasound frames recorded consecutively during a scan: we refer to them as the pre- and post-deformation frames. A window is placed around a point of interest in the pre- deformation frame, and the closest match in the post-deformation frame is located, hi practice, this is an optimisation problem, where the peak must be found in some suitable measure of signal similarity, such as the correlation coefficient (M. A. Lubinski, S. Y. Emelianov, and M. O'Donnell, Speckle tracking methods for ultrasonic elasticity imaging using short-time correlation, IEEE Transactions on Ultrasonics, Ferro- electrics, and Frequency Control, 46(1): 82-96, January 1999; J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li. Elastography: a quantitative method for imaging the elasticity of biological tissues. Ultrasonic Imaging, 13:111-134, 1991), sum of absolute (SAD) or squared (SSD) differences (S. Langeland, J. d'Hooge, H. Torp, B. Bijnens, and P. Suetens. Comparison of time-domain displacement estimators for two- dimensional RF tracking. Ultrasound in Medicine and Biology, 29(8): 1177-1186, 2003, R. L. Maurice and M. Bertrand. Lagrangian speckle model and tissue-motion estimation theory. IEEE Transactions on Medical Imaging, 18(7):593-603, July 1999, F. Viola and W. F. Walker. A comparison of the performance of time-delay estimators in medical ultrasound. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 50(4): 392-401, April 2003) or mutual information (M. I. Miga. A new approach to elastography using mutual information and finite elements. Physics in Medicine and Biology, 48(l):467-480, January 2003). Numerous phase-based approaches have also been developed (X. Chen, M. J. Zohdy, S. Y. Emelianov, and M. O'Donnell. Lateral speckle tracking using synthetic lateral phase. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 51(5):540-550, May 2004; M. O'Donnell, A. R. Skovoroda, B. M. Shapo, and S. Y. Emelianov. Internal displacement and strain imaging using ultrasonic speckle tracking. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 41 :314-325, May 1994; A. Pesavento, C. Perrey, M. Krueger, and H. Ermert. A time efficient and accurate strain estimation concept for ultrasonic elastography using iterative phase zero estimation. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 46(5): 1057-1067, September 1999), which exploit a property of the cross-correlation function peak, and are advantageous because of relatively low computational cost. Whichever technique has been used to match the windows, it is assumed thereafter that the mechanical displacement of tissue at the centre of the window is equal to the optimal window displacement. Window-matching is applied at positions throughout a grid over the acquired frame of ultrasound data, constructing a fine map of the displacement field.

A strain image can be produced by displaying spatial derivatives from the estimated displacement field. We consider in detail the problem of axial strain image foπnation, but the principles we derive are generally applicable. Strain estimation may be regarded as a stochastic process in which case the terms "mean squared error", "estimation noise" and "estimation variance" may be used interchangeably when referring to the typical discrepancies between actual deformations and the estimates that are recorded and displayed. Errors in strain images arise mostly from two sources. The first is displacement estimation error, which is well understood. Following Carter (G. C. Carter, Coherence and time delay estimation. Proceedings of the IEEE, 75(2):236-255, 1987) and (W. F. Walker and G. E. Trahey, A fundamental limit on the performance of correlation based phase correction and flow estimation techniques. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 41(5): 644-654, September 1994) has become popular to apply Cramer-Rao lower bound analyses (and variations thereon) to signals with known properties, thereby identifying a lower bound on the displacement estimation variance that could be achieved by a maximum likelihood estimator (F. Viola and W. F. Walker, A comparison of the performance of time-delay estimators in medical ultrasound. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 50(4): 392-401, April 2003).

The second source of noise has, however, been overlooked. The problem is estimation location variance: it is not generally true that the displacement estimate most closely tracks the actual displacement at the window centre. It was noted in an earlier study by (I. Cespedes and J. Ophir, Reduction of image noise in elastography. Ultrasonic Imaging, 15:89-102, 1993) that if there is intra-window compression and the signal envelope is not constant, then the actual estimation location is skewed towards higher amplitude portions of the widowed signal. This causes artefacts at the boundaries between regions of differing echogenicity, as demonstrated by Figure 1.

Figure Ia shows a B-mode image of RF data from a scan of a human arm. The signal is temporally compressed to simulate a uniform compressive strain of 1%. On a linear scale from black (0% strain) to white (2%), this should produce a uniform strain image with extremely low estimation noise, since the signal SNR is higher than could possibly be achieved in a real compression scan. However, Figure Ib shows that the standard correlation coefficient maximiser produces a strain image that is severely degraded (and misleading) owing to the AM (amplitude modulation) effect (which we explain later), while Figure Ic shows the (near perfect) result from applying the correction technique we describe later (strain estimation for both images used windows of length 13.5λ).

It is helpful at this point to introduce some of the terminology generally used in ultrasound imaging. An ultrasound imaging system generally employs a one- dimensional or two-dimensional ultrasonic transducer array (although sometimes only a single transducer may be employed), the array comprising typically 20 to 256 transducers in each dimension. Each transducer acts as both a transmitter and a receiver. The transducers are generally driven by a pulse of RF energy, typically in the range 1-20 MHz; the signal may be considered narrow band in the sense that a pulse is sufficiently long to include a number of RF wavelengths thus having a relatively well- defined frequency. The ultrasound transducer array is usually coupled to the tissue under investigation by an ultrasound gel or water; typically the ultrasound penetrates a few centimetres, for example up to 25cm, into the tissue under investigation, and the transducer array scans a region of a few centimetres in a lateral direction. The axial resolution is generally much greater than the lateral resolution, for example of the order of 1000 samples (in time) as compared with of the order of 100 lines laterally. So-called A-lines run actually from each transducer into the tissue under investigation; a so-called B-scan or B-mode image comprises a plane including a plurality of A-lines, thus defining a vertical cross section through the tissue. A two-dimensional transducer array may be used to capture perpendicular B-scan images, for example to provide data for a three-dimensional volume. A captured image is generally built-up by successively capturing data from along each of the A-lines in turn, that is by capturing a column of data centred on each ultrasonic transceiver in turn (although beam steering may be employed). However, when capturing data from a particular line, preferably a set of the transducers is driven, with gradually increasing phase away from the line on which the transducer is centred so as to create an approximately spherical ultrasonic wavefront converging on a focus on the line under investigation. The signals received from the transducers are summed with appropriate amplitude and phases to reconstruct the line data. This provides an RF (radio frequency) output which is usually time-gain compensated (because the amplitude of the received signal decreases with increasing probed depth) before being demodulated, optionally log-weighted and displayed as B-scan. Often the RF data is digitised at some point in the processing chain, for example prior to the demodulation, the remainder of the processing taking place in the digital domain.

In embodiments of the techniques we describe later both signal amplitude and signal phase are employed and therefore preferred examples of an ultrasound imaging apparatus suitable for implementing embodiments of the invention employ a pair of analogue-to-digital converters to provide in-phase and quadrature digitised signal components so that phase data is available.

An outline block diagram of an ultrasonic imaging system suitable for implementing embodiments of the invention is shown in Figure Id. We will describe how at least one-dimensional image data captured by a pulse-echo technique, in particular an ultrasonic imaging system, can be processed to determine improved strain field data, as shown in Figure Ic. The ultrasonic image data which is processed in embodiments of the technique comprises the digitised RF signal data, as shown in Figure Id, optionally with pre-processing in the analogue domain. (Where such pre-processing is employed it may take many forms, one example of which is illustrated). Broadly speaking the demodulated data may be processed by envelope detection and log weighting to provide a B-mode display and/or strain determination may be employed to provide a strain display. The demodulation illustrated in Figure Id extracts the amplitude (envelope) and phase information of the RF signal in a conventional manner and in some preferred embodiments the signal is digitised after demodulation so that there is no need to run the A/D at (twice) the RF frequencies employed - and thus in these embodiments the processed RF signal comprises a demodulated base band signal, hi other systems the RF signal may be digitised prior to demodulation.

A digitised I and Q (in-phase and quadrature) signal is frequently available in conventional ultrasonic imaging equipment and, conveniently, embodiments of the invention may be implemented by processing this signal using a suitably programmed general purpose computer or digital signal processor (DSP) and/or by using dedicated hardware.

In the description which follows it is often convenient to refer to samples in time rather than explicitly to position data, but typically a direct relationship is assumed between these two variables (position = velocity x time) assuming a typical speed of sound. For human tissue the speed of sound is normally taken as 1540 m/s (the speed for water at 49°C), which is accepted as a good estimate; for other materials other speeds may be employed. Similarly in the discussion which follows we will generally refer to axial strain (because the resolution of a system is typically highest in the axial or A-line direction) but it will be appreciated that the techniques we describe are also applicable to lateral strain in one or two dimensions (with a two-dimensional array), albeit normally with reduced precision because of the reduced number of samples).

We now continue, first with a discussion of the origin of the degradations which can be seen in Figure Ib. Broadly speaking this arises from the assumption that the position associated with the estimated displacement is the centre of a window. However, close inspection of Figure Ia reveals, for example, a diagonal stripe 10a which indicates a region of locally relatively high RF signal strength, and this corresponds with a similar diagonal strip 10b in the strain image of Figure Ib. As will be explained in more detail later, the stripe in Figure Ib in effect results from an incorrect association of the displacement estimation location with the centre of a window. A better estimation might be, for example, the highlighted region of the stripe. This example illustrates how degradation can arise. We now consider more closely some of the underlying theory.

It is observed that strain estimates are corrupted by an amplitude modulation (AM) effect, hi fact, the AM effect also degrades strain estimates within regions that are isoechoic, since the signal returned from a fine scatterer distribution dopes not have a constant envelope. Nevertheless, the AM effect is most serious in anisoechoic regions, where AM noise correlates strongly with the features in B-mode images, and can easily lead to severe misinterpretations of strain images.

It will be shown in the following section that the AM effect is often the primary source of error in ultrasonic strain images where it is not corrected. Two correction techniques were proposed by Cespedes and Ophir [ibid]. Firstly, log compression of the signal envelope reduces amplitude fluctuations, thereby shifting estimation locations towards the window centres. This is an effective means of mitigating the AM effect, and has consequently been reapplied in more recent strain imaging systems. The second suggestion was adaptive stretching, which compensates for intra-window compression by stretching the signal to enable a close match to the true displacement at all points. This has been shown to be a good way of reducing strain estimation noise, although such techniques are computationally expensive. The estimation location variance can also be reduced by using shorter estimation windows, but this is inevitably accompanied by reduced accuracy in the displacement estimates, since displacement estimation variance increases as the reciprocal of the window length.

The AM effect is present in all displacement tracking methods that use amplitude information, including methods based on the (normalised) correlation coefficient. To eliminate the AM effect, the amplitude must be entirely suppressed, as in one-bit compression, but this can bring unwanted side effects. We examine the AM effect from a theoretical standpoint, leading to a surprisingly simple AM correction method (AMC). Experiments have been performed using simulated RF ultrasound data to compare the performance of phase and correlation coefficient methods, and to evaluate the efficacy of correction by AMC, log compression and one-bit (limiting) compression in both cases. All of the corrections are computationally efficient and suitable for use in real- time imaging systems. Further experiments were performed using a direct strain estimator with adaptive stretching, which is slower but provides an AM suppression benchmark, which in fact embodiments of our techniques better.

Amplitude modulation theory

We analyse the estimation of strain from a set of window displacement estimates. For the sake of clarity, we examine the simplest method for converting ID displacement estimates to ID strains, by taking the difference between displacements at consecutive windows, and dividing this by the spacing between the assumed estimation locations.

S = ^Zώ. (1)

T2 ^~ X\ s is the strain estimate, J₁ and J₂ are the displacement estimates for windows 1 and 2 respectively, and f_x and f ₂ are assumed to be the estimation locations. It is commonly assumed that Equation 1 contains only two random variables: J₁ and J₂. Here we examine the neglected variables, f ₂ and f , . New variables D and F are defined to simplify the strain calculation.

D = d₂ -l (2)

F = -±— (3)

s = DF (4)

The sources of estimation noise are illustrated in Figure 2. We will assume that errors in D and in F are uncorrelated. This allows the strain estimation variance,

, to be expressed in a simple form (in which the first and third terms comprise AM).

μ^ is the expectation of D , which for an unbiased estimator is equal to the actual difference, D , between the displacements of the two windows, σ^ is the variance of D , which is approximately equal to the sum of the variances of the individual displacement estimates, J₁ and J₂ (it is exactly equal only if errors in J₁ and J₂ are uncorrelated, which is not the case for overlapping windows), μ^ is the expectation of the reciprocal location spacing estimate, F , which may correspond to the reciprocal of the spacing between consecutive windows. Finally, σ| is the mean squared error between F and the actual reciprocal spacing, F . In general, F is not equal to the reciprocal of the window spacing, since the actual estimation locations, τ₂ and τ_λ , do not generally correspond to the window centres.

We want to know what impact the terms in Equation 5 have on strain image quality. We consider a quality measure denoted SNR₆ , which has previously been defined (I.

Cespedes and J. Ophir. Reduction of image noise in elastography. Ultrasonic Imaging, 15:89-102, 1993; T. Varghese and J. Ophir. A theoretical framework for performance characterization of elastography: the strain filter. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 44(1): 164-172, January 1997) and can be measured experimentally in images where the underlying strain field is known to be homogeneous.

μ_s is the mean strain estimate and σ- is the standard deviation. The presence of μ². in the third term of Equation 5 becomes important when SNR₆ is evaluated. The noise contribution from the AM effect is therefore proportional to the strain, s , so the AM effect is expected to become the dominant source of strain estimation noise as the level of strain increases.

Equation 7 is derived by substituting the RHS of Equation 5 into Equation 6. The final result incorporates some simplifying assumptions. (1) μ. = s . (2) The assumed value of

F is usually a constant, i.e. μ- = F = K_x . (3) μ_ά = K₂S where K₂ is a constant (the expected shift equals the strain multiplied by the window spacing). We next describe some examples with pulse train signals.

Window matching tracks the displacement of the enclosed signal. However, if displacement varies within the window, then the actual signal displacement cannot be matched at all points. The location at which the actual displacement of the signal is equal to the displacement estimate varies depending on both signal and displacement field properties. In general, the estimation location comes from a random distribution throughout the window. It has low probability density at the ends, and in the absence of additional information its expectation is the window centre. Where the location cannot be estimated, it is best to assume that windows sample displacement at their centres. Unfortunately this means that the AM effect introduces displacement and strain estimation noise, as illustrated in Figures 1 and 2.

It is not possible to devise an estimator that both produces optimal displacement estimates and samples displacement at the centre of the window. This is because some portions of the signal may contain no information, or the quality of the information may be variable. This is demonstrated by examples with pulse train signals in Figure 3, which shows extreme examples of the AM effect, showing the output of a perfect displacement estimator operating on different pulse train signals with uniform strain. The strain (displacement gradient) is underestimated in Figure 3 a and overestimated in Figure 3b. In the absence of information between the pulses, an optimal displacement estimator tracks the displacement of the pulse(s) within each window. The example medium has been deformed by a uniform strain field, so displacement varies linearly with distance. The assumption of estimation at the window centre now leads to significantly different strain estimates if (a) overlapping windows track the same pulse, or (b) neighbouring windows track pulses at their extremities. When a uniform strain, s , is being tracked, and there is no displacement estimation error, the AM effect nonetheless distorts the result, such that the strain estimation lower bound is 0 for overlapping windows, and the upper bound is s x ¹^- . T is the window length and At is the spacing between successive windows. For non-overlapping windows the lower bound is s x ^f- . Referring again to Figure 3, in effect only a short part of the window has substantially all the displacement data. Since strain is displacement divided by the difference in locations the point at which the displacement is measured matters. Although Figure 3 is exaggerated this type of effect is particularly pronounced at, for example, the boundary between fat and muscle or other tissue types.

A real ultrasound signal is not generally a pulse train or the AM effect could be corrected by noting that displacement estimation occurs at the pulse locations. However, real ultrasound signals do incorporate amplitude variations, which are often large even over small distances. Lower amplitude sections usually have lower SNR, and a good displacement estimator should incorporate a mechanism for preferentially tracking the most reliable data. Ideally it should also be possible to estimate the actual displacement location when this is not equal to the window centre.

Summary of the Invention

According to a first aspect of the invention there is therefore provided a method of processing at least one-dimensional image data captured by an imaging technique to determine deformation data defining an at least one-dimensional deformation in an imaged object, said deformation data comprising at least one data pair, said data pair comprising displacement data and corresponding displacement location data, the method comprising: inputting first and second sets of said image data, a said set of image data comprising imaging signal data including at least one of signal magnitude data and signal phase data for an imaged region of said object; positioning a first window at a first position on said first set of image data and determining a corresponding position for said first window in said second set of image data; determining first displacement data defining a first displacement estimate from said positions of said first window in said first and second image data; and determining an estimated location of said first displacement estimate responsive to said imaging signal data within said first window for at least one of said first and second sets of image data; whereby said at least one data pair comprises said first displacement data and said estimated location of said first displacement estimate. The deformation data may define, for example a displacement or strain field within the object (the sets of image data will then generally correspond to different deformations of the imaged object). However the deformation data need not be employed to determine a strain estimate per se as there are many ways in which this data may be processed (or it may simply be displayed in a substantially unprocessed form).

In some preferred embodiments the technique is applied to a plurality of windows, for example at successive positions in the one (or more) dimensional image. The window positions may overlap, as described further later. In some applications, however, a single window position may suffice. For example where tissue is imaged there may be zero motion at the probe surface, which may be taken as a reference to provide, say, an estimate of a mean strain between the probe surface and the window position.

In some preferred embodiments therefore the deformation data comprises a set of data pairs, and the method further comprises positioning a second window at a second position on said first set of image data and determining a corresponding position for said second window in said second set of image data; determining second displacement data defining a second displacement estimate from said positions of said second window in said first and second image data; and determining an estimated location of said second displacement estimate responsive to said imaging signal data within said second window for at least one of said first and second sets of image data. Then a second data pair of the set of data pairs may comprise the second displacement data and the estimated location of said second displacement estimate.

Preferably the image data comprises data captured by a pulse-echo imaging technique such as ultrasonic imaging. However embodiments of the technique may also be applied to CT (computer tomography) elasticity imaging and even, for example, to optical imaging techniques, for example for inspecting strain in skin.

In preferred embodiments the at least one-dimensional image data preferably comprises digitised RF (radio frequency) data, either before or after demodulation. In some preferred embodiments this data is in the form of in-phase and quadrature digital signals, although other data formats may also be employed. Typically one of the sets of image data defines a pre-deformation frame (here frame including one-dimensional data) and the other a post-deformation frame. It will be understood that either of these may be considered as a reference (depending upon whether positive or negative strain is considered) for the estimated displacement estimate location.

Once the deformation data has been obtained it may be used in any convenient manner; there are many ways in which such data may be used. Typically information derived from this data is displayed to an operator of the system, for example as an image of a displacement or strain field in the imaged object. Generally the deformation data defines strain at a plurality of locations within the imaged object although, at its simplest, only a single strain value for the imaged object may be defined.

Preferably a displacement estimate location for displacement data determined from a window is determined dependent upon the signal magnitude (envelope) data for corresponding windows in both the first and second sets of image data, in embodiments by determining a centroid using the signal envelopes in corresponding windows. The centroid may be determined by evaluating SUM {weightings x locations} / SUM {weightings} . More particularly the method determines the centroid of signal magnitude squared when using only one set of image data or frame for the estimate or, in some preferred embodiments, the centroid of the product of the pre- and post- deformation signal magnitudes, that is using the signal magnitude data for both the first and second sets of image data or frames. (Some benefit may be obtained by using the signal magnitude data from just one of a pair of corresponding windows). The determining of an estimated displacement estimate location may additionally (or alternatively) be responsive to signal phase, more particularly a relative phase difference between signals for two corresponding windows (a difference in signal phase between the first and second sets of image data).

Broadly speaking in embodiments the determining of a displacement estimate location involves establishing an approximation of the following form: estimate for d = (SUM: weightings*displacements) / (SUM: weightings), where the sum is over a window. Optionally this may also take into account signal phase. Having established an approximation of this general form, (time or position) data can be included in the numerator and this leads to an expression for the determination of a displacement estimate location as shown, for example, in equations (28) and (31) later. In general, any signal property that can be measured may be incorporated in such a weighting approximation; signal magnitude may just be one component of the weighting expression. For example, for correlation function methods such as EPZS (Efficient Phase Zero Search) one technique is to use the envelope power, weighting = envelope in first window * envelope in second window, but a more accurate approximation is afforded by including phase, so that, for example: weighting = envelope in first window * envelope in second window * sin(phi)/phi where phi = phase in second window - phase in first window - phase of overall correlation function, hi embodiments, in particular where the correlation function has zero phase, the value of phi may be determined by a difference in phase of the signals at corresponding points in the two windows. In more mathematical terms, the weighting, W, may be given by E₁ (m) ■ E_x (m) ^■ sin(Φ(m)) / Φ(m) where Φ = (phase of Ei at point (sample) m - phase of

E₂ at point (sample) m) and where the phase of Ei and E₂ may be taken directly from the digitized I and Q signals.

An example of an AM correction technique in which only phase information may be used is in the EPZS_L2 (Efficient Phase Zero Search, limiting log compression) procedures, described in more detail later, in which the magnitude information has been discarded. In such a procedure the weighting, W = sin(Φ)/ Φ may be employed for measurable performance improvements. Another heuristic example of phase weighting is ((7T - Φ)/7r)ⁿ, which deweights large phase differences (large Φ) at which, for example, phase wrapping errors are more likely. The degree of phase (de)weighting may be detrmined by selecting n, for example n = 0 corresponding to no phase deweighting and n = 3 corresponding to severe phase deweighting.

Thus in embodiments phase data alone may be used to determine a displacement estimate location. In embodiments the determining of the window positions may include interpolating between positions corresponding to (digital) data samples. For example using the gradient at the zero phase point a window position may be determined to an accuracy of better than Vi₀, V₁₀₀ or Viooo of a sample (for example, dependent on the level of decorrelation). In effect locating the matching window involves linearly interpolating between samples, for example envelope values of samples, to determine a more accurate position. The techniques we describe generally involve cross-correlating the received image data to determine corresponding window positions in the two sets of data. They may be applied to cross-correlation methods which determine the maximum value of a cross-correlation or, more preferably, of a cross-correlation coefficient, and to phase- based methods, which estimate a window position from the phase of a value of the cross-correlation function. However there are several other techniques which may also be employed to locate a matching window, as mentioned in the introduction.

Generally successive pairs of windows (that is pairs comprising corresponding windows in the two sets of image data) are processed, moving down the A-lines (and optionally back up the A-lines to help reduce phase ambiguities) away from an ultrasonic transducer or transducer array. In embodiments the captured data comprises data from a two-dimensional transducer array and two-dimensional deformation data is determined.

In embodiments the deformation data may be used to determine an actual strain for the object or to infer or image a property of the object such as elasticity (including one or more viscoelastic moduli). However in other applications a greyscale or colour image of the raw data may simply be displayed. The embodiments of the technique we describe are sufficiently sensitive to rely upon stresses induced by an operator's hand to generate a strain field. However a device such as a controlled vibration device may be employed to provide a controlled stress to the imaged object and this stress data may then be combined with deformation, in particular strain field data determined using the above described techniques in order to calculate and/or display/image one or more properties of the imaged object, such as elasticity. There are other techniques for stressing the imaged object, for example using a beam of (focussed) ultrasound to induce an internal mechanical stress. The stress from physiological motion may also be employed, for example, stress due to blood pressure variations during the cardiac cycle. In preferred embodiments the object comprises biological tissue, preferably living human or animal tissue although other biological material such as foodstuffs, for examples meat or fruit may be imaged. Although preferred embodiments of the method are employed in the field of ultrasonic imaging, applications of the technique are not limited to this field, hi particular the method may also be applied to magnetic resonance imaging (MRI) in which case the pulse comprises an RF electromagnetic pulse and the echo a spin-echo. It is known to employ MRI for elasticity imaging (sometimes referred to as MRE - magnetic resonance elastography) - see for example, Oida, T., Amano, A. and Matsuda, T, "MRE: In vivo measurements of elasticity for human tissue", Informatics Research for Development of Knowledge Society Infrastructure Conference, 1-2 March 2004, p.57-64 - and the techniques we describe may also, in embodiments, be applied to MRE. hi still further embodiments the techniques may be applied to CT elasticity imaging.

The invention further provides processor control code to implement the above-described methods, for example on a general purpose computer system or on a digital signal processor (DSP). The code may be provided on a carrier such as a disk, CD- or DVD- ROM, programmed memory such as read-only memory (Firmware), or on a data carrier such as an optical or electrical signal carrier. Code (and/or data) to implement embodiments of the invention may comprise source, object or executable code in a conventional programming language (interpreted or compiled) such as C, or assembly code, code for setting up or controlling an ASIC (Application Specific Integrated Circuit) or FPGA (Field Programmable Gate Array), or code for a hardware description language such as Verilog (Trade Mark) or VHDL (Very high speed integrated circuit Hardware Description Language). As the skilled person will appreciate such code and/or data may be distributed between a plurality of coupled components in communication with one another.

In another aspect the invention provides apparatus for processing at least one- dimensional image data captured by an imaging technique to determine deformation data defining an at least one-dimensional deformation in an imaged object, said deformation data comprising at least one data pair, said data pair comprising displacement data and corresponding displacement location data, the apparatus comprising: an input for first and second sets of said image data, a said set of image data comprising imaging signal data including at least one of signal magnitude data and signal phase data for an imaged region of said object; a system for positioning a first window at a first position on said first set of image data and determining a corresponding position for said first window in said second set of image data; a system for determining first displacement data defining a first displacement estimate from said positions of said first window in said first and second image data; and a system for determining an estimated location of said first displacement estimate responsive to said imaging signal data within said first window for at least one of said first and second sets of image data; and wherein said at least one data pair comprises said first displacement data and said estimated location of said first displacement estimate.

These and other aspects of the invention will now be further described, by way of example only, with reference to the accompanying figures in which:

Figures Ia to Id show, respectively, a B-mode image of a human arm, a conventionally generated strain image, a strain image generated using an embodiment of a method according to the invention illustrating a near-perfect result, and a block diagram of an ultrasonic imaging system for implementing embodiments of the invention;

Figure 2 shows sources of error in strain data estimation from captured RF image data;

Figures 3 a and 3b illustrate the origin of an amplitude modulation (AM) effect;

Figure 4 shows an example RF ultrasound signal;

Figures 5a and 5b show, respectively, a procedure for determining estimated displacement and displacement location data from ultrasound scan data, and a procedure for displacement location correction according to an embodiment of the invention;

Figure 6 shows an example B-scan image from simulated RF data; Figure 7 illustrates displacement offsets in adaptive strain estimators;

Figure 8 shows signal-to-noise ratio against window length for Efficient Phase Zero Search (EPZS)-based techniques;

Figures 9a and 9b show, respectively, signal-to-noise ratio against window length for CCM (Correlation Coefficient Maximiser) and EPZS and ASE (Adaptive Strain Estimator)-based techniques;

Figures 10 to 1Od show strain images using, respectively EPZS (SNR_e = 1.63), CCM (SNR₆ - 1.62), EPZS_A (SNR_e = 3.86), and CCM_A (SNR₆ = 2.05);

Figure 11 shows SNR₆ against compression factor for EPZS and CCM-based techniques at 0.01% strain;

Figure 12 shows SNR₆ against compression factor for EPZS and CCM-based techniques at 0.5% strain;

Figure 13 shows SNR_e against compression factor for EPZS and CCM-based techniques at 4% strain;

Figure 14 shows SNR_e-strain characteristics for EPZS-based techniques;

Figure 15 shows SNR_e-strain characteristics for CCM-based techniques;

Figure 16 compares SNR_e-strain characteristics for EPZS, CCM and ASE-based techniques;

Figure 17a-e show an imaged olive/gelatine phantom showing (a) B-scan, (b) EPZS (white=255=l% strain), (c) EPZS-Ll, (d) EPZSJL2, (e) EPZS-A; Figures 18a-e show a gelatine phantom with two regions, showing (a) B-scan, (b) EPZS (255=0.8%), (c) EPZS_L1, (d) EPZS JL2, (e) EPZS_A;

Figures 19a-e show human male breast tissue showing (a) B-scan, (b) EPZS (255=0.8%, (c) EPZSJLl, (d) EPZS_L2, (e) EPZS_A; and

Figures 20a-e show human male calf images showing ((a) B-scan, (b) EPZS (255=0.8%, (c) EPZS_L1, (d) EPZSJL2, (e) EPZS_A.

Broadly speaking we will describe a technique to approximate an estimator as a weighted sum of displacements throughout a window, and use the spatial distribution of those weightings to determine a location where the estimate is valid. This can be done by taking the centroid of those weightings, in effect to estimate the location at which the displacement estimate is valid.

We first describe an analytical investigation of the AM effect with particular reference to phase-based methods.

We derive an approximate expression for the AM effect when windows are matched by identifying the zero crossing of the complex cross-correlation phase. Phase-based methods operate on analytic signals with real and imaginary parts, which are produced by applying the Hubert transform (or some approximation thereof). The complex cross- correlation function and its phase may be expressed as follows.

(_aι,a₂)(nAt,d) = "γ a, (t)a₂ (t + d) (8) t=nhl

Ci_x and α₂ are analytic ultrasound signals, ^* denotes the complex conjugate, ntύ is the location of the beginning of the analysis window in the pre-deformation signal, T is the window length, and d is the candidate displacement applied to the post-deformation window to look for a match. Eventually the match or displacement estimate, J_n , is found where φ has a zero crossing. Broadly speaking the technique is to slide the window along until a zero-crossing is detected. Preferably (as described, for example, in US 6,520,913) when looking for the window match in the post-deformation frame the sliding of each successive window starts from a position which is determined by the pre-deformation window position plus the sum of the displacements of previous post- deformation windows along the A-line (so that the searching does not get progressively longer).

It will be noted that if φ is only expressed in the range [-π, +π] then a zero crossing occurs on average once for every wavelength shift in d . hi some implementations it may therefore be useful to incorporate a system for guiding the search to help ensure that the correct zero crossing is selected, for example by starting at the top and at the bottom of a column, and then taking the better of the two results. Without this an error in matching a window early in a column (A-line) can cause the rest of the matching to fail, showing a line error in the image - "dropout" (visually untidy but not directly relevant to the techniques we describe here). This is analogous to eliminating "peak- hopping" errors from correlation coefficient analysis and there are many conventional techniques for addressing this (for example, Walker & Trahey, ibid). Other error detection and correction techniques are described in (J. E. Lindop, G. M. Treece, A. H. Gcc, and R. W. Prager, 3D elastography using freehand ultrasound, 2006, Ultrasound in Medicine and Biology; Y. Zhu and T. J. Hall, A modified block matching method for real-time freehand strain imaging, Ultrasonic Imaging, 24:161-176, 2002).

To analyse the properties of phase-based methods, we use a simple signal model with no noise, where decorrelation occurs only as a result of the ID signal stretching that accompanies mechanical strain. Our model of the pre-deformation signal, a_λ , is a constant frequency sinusoid, scaled by a positive real signal envelope, / . This is illustrated in Figure 4, which shows a signal model comprising a constant frequency sinusoid with an arbitrary signal envelope, subject to an arbitrary deformation. Real ultrasound signals are not dissimilar to Figure 4. a_x{t) = f(t)e^jωl (11) The main limitation of this model is the constant frequency assumption, but real RF ultrasound signals are narrowband, although the frequency may be substantially constant over short distances.

hi our model the post-deformation signal, a₂ , is produced by an arbitrary temporal warping of a_x , such that every point, α, (t) , undergoes a displacement, d{t) . a₂(t + d(t)) = a₁(t) (12)

This is a simplification of the signal transformation that occurs in a real compression scan. Firstly, it will be noted that a uniform strain in our model gives rise to a change in the frequency centroid in the post-deformation signal, which will not usually be reflected in reality (although local frequency changes do occur). Secondly, we have assumed that the only change to the signal envelope will be a ID warping. In reality, changes in the interference patterns of closely spaced scatterers introduce unpredictable components in the post-deformation signal, resembling the addition of an uncorrelated narrowband noise signal. Furthermore, axial compression in real materials with finite compressibility is inevitably accompanied by additional motions in the lateral and elevational directions. Nevertheless, we pursue analysis with our simplified model, and the predictions, later tested on real and simulated ultrasound data, produce good results.

We examine the properties of the signals in matched windows. In general, the estimated displacement is similar but not equal to the local displacement at each position in the window. We therefore introduce a new symbol, t₂ , denoting the pre-deformation location in α, , of the data with which Ci₁ (t) is matched. t₂ + d(t₂) = t + d_π (13)

The complex cross-correlation function at the match is now expressed as follows. ιiAt+T

{a_va₂)(nM,l) = ∑ a_λ {t) a₂ (t + l) (14) t-nAt nΔt+T

= ∑ <(0 ^α. ('₂) (15)

I=IiAt

m/(Φ^Mh-^ή (i6)

In order to satisfy the match criterion (Equation 10), the imaginary part of the complex cross-correlation function must be zero. nAt+T

3 ∑ f(β)f (t₂)e .Mh-ή = 0 (17) t=ιιAl

This leads to an alternative expression for the phase zero condition. nAt+T

∑ f(t)f {t₂)sm(ω(t₂ -t)) = 0 (18)

I=IiAt

It is noted that t₂ - 1 = ^_n - d(t₂) is the local discrepancy between the displacement estimate and its actual value. This is small, so ω(t₂ -t) «: -f- at all points in the window for typical window lengths and operating strains. It follows that Equation 18 may be simplified by applying the small angle approximation. iiAl+T

∑ f(t)f (t₂)ω(t₂ -t) = O (19)

I=ItAt

Equation 19 can be converted to an expression with clearer relevance to the physical deformation by examining the term t₂ - 1 . This is performed as follows, employing the relation from Equation 13, and expanding a Maclaurin series about d(t) . t₂ -t = Q_n - d(t)} - {d {t₂ )- d(t)} (20)

= Q_a -d(t)} -^P- {t₂ -t} - O{(t₂ -t)²} (21) dt

= {l -d(t)}-_S{l -d(t₂)}-θ{(l -d(t₂)f} (22)

In developing this theory second order terms will be neglected, as will the term scaled by s (strain), since the vast majority of previously documented ultrasonic strain imaging systems operate with s « 1.0 (although the results are not limited by this). Typical strains are < 10% and strains of, for example, 0.05% induced by vibrations of an operator's hand, can produce useful images. Now the result from Equation 22 is substituted into Equation 19. nAt+T

∑ f(t)f [t₂)ω(ά_n - d(t)) = ^Q (23)

Rearrangement yields a good approximate formula for the displacement estimate, _j .

We have shown that an approximation of the phase-based displacement estimate is a weighting of point displacements by the cross power of the local signal envelope.

Amplitude modulation correction

We now show that the actual estimation location can be estimated for the important case where strain may be considered constant at the scale of the individual windows (although we have found, in practice, that the results we describe are more generally applicable). The constant strain condition is expressed mathematically as follows. d(t) = a + st (25)

We substitute this into Equation 24, and rearrange to produce a convenient form for the approximation.

The location estimate, f „ , is defined to be the position (or equivalently time t) at which the displacement estimate approximation ( ^_n ) is equal to the actual displacement (d), i.e. ^_Λ = α + Jf „ . Hence,

In equation (28)/(^ and /^ are the envelopes of the received ultrasound signal in the pre- and post-deformation frames, which may be obtained from the magnitude of the digitized 1 and Q signals, either before or after demodulation. The sum is over the respective windows and, in effect, calculates the centroid of the product of these two signals. In embodiments, however some benefit may nonetheless be derived from using a simplified approximation to equation (24), for example using only one or the other envelope.

The location estimates of equation (28) are substituted into Equation 1 to refine the strain estimates. In effect, under the assumption of constant strain equation (28) gives the position at which the displacement estimate is not accurate. Broadly speaking, the displacement estimate is weighted at each location within a window by a combination of the RF signal envelopes in a window and its corresponding window on the pre- and post-deformation frames. This amplitude modulation correction (AMC) also allows a more accurate identification of the image region corresponding to the space between successive displacement estimates, thereby producing a more accurate correspondence between the physical locations of tissue features, and their apparent locations in strain or displacement images.

More complex corrections may be employed. For example in a case where successive windows overlap (although not necessarily limited to this case) equation (25) may be replaced with a higher order polynomial fit to the deformation field. Then, solving the simultaneous equations this produces (in place of equation (26)) can potentially be used to produce a more accurate estimate. Such a technique is potentially useful in a displacement field where higher order derivatives are highly significant.

Thus in a linear approximation of displacement (equation 25 above) we effectively determine an intercept of the actual displacement on this line. However one or more additional terms may be included in a polynomial approximation along the lines indicated below:

∑ W(f)t ∑ W(t)t² d - CC_n + a, ^■ + a. + ...

Σ W if) ∑ W(ή

In such a case the location of the displacement estimate may be determined by the intercept of the actual displacement on this curve. However, equivalently in embodiments of the invention one or more of the parameters of the above expansion may be used directly (or indirectly) to provide the deformation data. For example, the term a_λ in the above expansion provides strain data directly. Thus the location of the displacement estimate may be considered as effectively comprising one or more of these parameters.

The above described approach can also be applied to correlation coefficient methods which have to date been the most popular approach for displacement tracking, at least for ultrasonic strain imaging. The correlation coefficient for real RF signals r_λ and r₂ at window n with a candidate shift d is evaluated as follows:

The displacement estimate is chosen to maximise the correlation coefficient. d_H = arg max /?,.,. (n&t, d) (30)

In common with the analysis of phase-based methods, it is highly desirable to derive a similar estimation location expression. The starting point is to identify the properties of stationary points (including the maximum) by differentiating p with respect to d .

Because it is difficult to derive an analytic expression for f „ in the case of correlation coefficient methods, we apply the following heuristic, which is motivated by an assumption that the AM effect on correlation coefficient methods is similar to the effect on phase-based methods, for which AMC has already been derived:

\ rΛt)rM + d ) \ t

∑_t,nAt \ Λ (t)r₂(t _{+ d)}) \ Simulation results are included later to show that this is a useful technique.

Referring next to Figure 5a, this shows a procedure which may be employed in embodiments of the invention for determining displacement estimation data and estimation location data (in time or equivalently, spatial position). Thus at step S500 receives scan data, that is RF signal data comprising magnitude and phase data, for two or more frames (in this context a frame of data may be one-dimensional, two- dimensional or even three-dimensional). Typically data for two B-scan images is input, pre- and post-compression, each comprising a plurality of A-lines (although at minimum a single A- line for each image will suffice). In a typical embodiment there may be 10 - 200 B-scans inputted per second and successive pairs of frames in time are considered as pre- and post-compression images (i.e. images 1 and 2, then images 2 and 3, and so forth). In embodiments a substantially real time strain data display may be provided.

Also at step S 500 the first window position is initialised on an A-line, generally beginning at a point in the imaged object adjacent the transducer. Preferably a window has a length of at least twice the wavelength of the ultrasound in the object. As indicated in step S502 the procedure steps a window in the first (pre-deformation, for example pre-compression) frame down the A-line, generally at regular intervals or according to a grid, and for each successive position looks for a match in the second (post-deformation, for example post-compression) frame; successive window positions may overlap. Thus, at step S 504 a window is positioned in the first frame and then, at step S 506, the procedure searches for a match for the window in the second frame, typically by cross-correlation. In preferred embodiments the procedure begins searching at an estimated displacement - for example for the second matching window position in the second frame, the procedure beings searching at the position in the first frame plus the displacement of the matching first window in the second frame plus the estimated displacement of the matching second window in the second frame. This type of iterative procedure is described in Pesavento (ibid) and US6,520,913 (hereby incorporated by reference). Once the matching window in the second frame has been found the procedure records the displacement (S 508), and then increments to the next window position (S510) the procedure then looping back to step S504 until there are no more window positions to process. At this point a set of window positions (time or spatial location) and corresponding displacement values has been determined. The position and displacement data may be processed in a number of ways or merely recorded, for example in a data file. Typically these data are employed to calculate strain values (S512) which are then converted to greyscale data and displays (S514), although sometimes the raw displacement data is displayed (in which case the image gets progressively darker down the screen). The calculation of strain values may comprise a straight forward application of equation (1) above, or may employ a more sophisticated technique, for example fitting a straight line to three, five, seven or more adjacent points, for example with a least squares fit, to determine a strain value. Other techniques may also be employed. The strain may be generated by small motions of an ultrasonic probe or a substantially known stress field may be generated automatically. There are many ways to generate such a stress field, for example controlled vibration of a probe, and with stress data elastic properties of an imaged object may be determined, at its simplest from a ratio of stress to strain. Other tasks may also be achieved, for example tracking of shear wavefronts for fast characterisation of tissue mechanical properties.

There are many variants of the above-described procedure to which embodiments of the technique described herein may be applied. For example, rather than positioning windows at regular intervals according to a grid window warping may be employed, for example using adaptive stretching Lagrangian estimation. The main point is, however, that there should be a set of estimated displacements and corresponding estimation locations to which embodiments of the technique described herein can be applied.

We next refer to Figure 5b, which shows a procedure for implementing an embodiment of a technique according to the invention. Although this is shown at a separate procedure, in practice it may often be convenient to run the steps of the procedure in parallel with steps S504 to S510 of the procedure of Figure 5a, that is in parallel with the window positioning and matching described above. Alternatively the procedure may be performed between steps S508 and S510 of Figure 5a, or between steps S510 and S512.

Referring to Figure 5b, inputs to the procedure comprise a set of displacement and displacement location data and either raw or demodulated RF signal data, preferably as in-phase quadrature digitised RF signal data, to provide the envelope and phase information shown in Figure 4 (step S550). The procedure then operates with successive pairs of windows in the first and second (pre-and post-deformation) frames mentioned with respect to Figure 5a above, generally employing at least the envelope of the RF signal data from each window of a pair (S 552). At step S554 the procedure calculates a weighted average estimated distance into the relevant window for the displacement estimate, that is a position within the window at which the displacement estimate is considered to apply, hi one embodiment the procedure determines

/V-I

* - ^■

More generally however the procedure may determine an estimated location for the displacement estimate based on a sum of a product of weightings and signal value samples, in particular based on:

/V-I

m - ^■ m=0

/V-I Σ W(m) m=0

hi the particular previous example, JV[TTl) — E [JTi)E [ Tn) . Some other example

weightings have been described above and employ the magnitude and/or phase information from one or both windows of a pair, hi general W(m) comes from an approximation for d as SUM (weightings * d) / SUM weightings.

In some preferred embodiments the weighting effectively comprises a centroid of the envelopes of both windows of a pair of matching windows in the pre- and post- compression frames. This is illustrated in step S554, where the sum is over the number of (A/D) samples, N in each window, the weighted average estimated distance into a window being calculated in terms of a number of samples (time, or equivalently spatial position). Optionally interpolation between samples may be employed for increased accuracy. Further, as previously mentioned, other weightings may be employed, for example in a simpler calculation the squared envelope of a single window of a pair. In the example step S554 the estimated distance through the window in terms of number of samples ( m ) defines a percentage of the linear length of the window (N samples) which is used to modify the relevant displacement location (S 556). The procedure then continues (S558) as before (steps S512, S514), making use of the corrected data in any desired manner, for example to display a corrected strain image of an object as shown in Figure Ic.

AMC increases the utility of displacement estimates from a spatially varying displacement field by estimating the actual estimation location. An alternative approach for handling the AM effect is to reduce the level of amplitude variation, for example by log compression of the signal envelope. This may be a useful technique in some circumstances, but the AM effect may actually be beneficial for high quality displacement estimation.

We describe later (under "Benefits of amplitude modulation correction") an analysis of a simple model of a generic displacement estimator, where short windows produce unreliable estimates, but the estimation variance can be reduced by using longer windows to take a weighted moving average. Following reasonable assumptions, it is shown that an optimal displacement estimator weights the importance of different signal sections in proportion with the local cross power, r, (t)r₂ (t + d_n) • This outcome is similar in form to the approximation in Equation 24 for phase-based methods. It implies that the weighting becomes suboptimal if the amplitude is compressed, thereby reducing the accuracy of the displacement estimator. We therefore expect that if location estimation such as AMC is performed accurately, then the lowest strain estimation noise is achieved in the absence of log compression. It is less clear how far these conclusions apply to correlation coefficient methods, but the correlation coefficient also incorporates a weighting of some form, since high amplitude sections within a window have a greater impact on the overall correlation coefficient value. We next briefly review adaptive strain estimators before discussing experimental methods and results.

Adaptive strain estimators work on the principle of reversing the deformation that has occurred, to obtain the best match to the pre-deformation signal. In ID, an adaptive strain estimator uniformly stretches the post-deformation window until its similarity to the pre-deformation window is maximised. If the local strain is actually uniform, adaptive strain estimation has the advantage of being able to correctly match the displacement at every point within the window. This means that for uniform strains the question of estimation location is ύτelevant, because the correct displacement can be found everywhere. Tests of adaptive strain estimation on uniform strain simulations are therefore expected to be independent of the AM effect, and it is for this reason that we employ an adaptive strain estimator as our AM suppression benchmark.

Experimental methods

Simulation

Simulated RF ultrasound data has been generated using Field II (J. A. Jensen. Field: a program for simulating ultrasound systems, In Proceedings of the lθ" Nordic-Baltic Conference on Biomedical Imaging, volume 4, pages 351-353, 1996). The simulations have 2x lO⁵ scatterers positioned at random according to a uniform distribution throughout a 50x 50x 6 mm volume, with random scattering strengths distributed uniformly over the range [0,γ_max] . The probe parameters model the 5-10 MHz probe of the Dynamic Imaging Ltd (UK) Diasus ultrasound machine, for which the point spread function has been measured experimentally — the pulse has a centre frequency of 6.0 MHz and bandwidth 2.1 MHz — and the sampling frequency is 66.7 MHz.

For each frame 128 A- lines have been simulated, spanning 40 mm in the lateral direction, recorded to a depth of 40mm. Simulations have been performed at a range of compressions (0%, 0.01%, 0.1% 0.5% 1.0%, 2.0%, 4.0%) by rescaling the axial spacing of the scatterers. This is important, because the relative performance of the strain estimation algorithms we compare is strain dependent. Five data sets have been generated for different scatterer fields. This contributes to the reliability of the results, which record the mean and standard deviation across the five data sets.

The Field II output was converted to the RF ultrasound format of the Stradwin freehand 3D ultrasound system (http://mi.eng.cam.ac.uk/^~rwp/stradwin/). RF samples are recorded with 16-bit signed integer precision. To ensure reproducibility of the resultant SNR and AM effects, the signals were normalised before conversion, such that in all cases the mean power is fixed at V_{1 m} = 2¹⁰ , corresponding to a mean SNR of 71 dB in the presence of quantisation noise. Tests have also been performed on simulated data with additive white Gaussian noise, reducing the SNR to 20 dB. Figure 6 shows an example B-scan from the simulated data.

In vitro and in vivo scanning

Scans have been performed using a Dynamic Imaging Diasus ultrasound machine with a 5-10 MHz probe, sampled at 66.7 MHz by a Gage CompuScope 14200 analogue-to- digital converter, with a PC running the Stradwin freehand 3D ultrasound software. As per previous work (Gage Applied Technologies Inc., Canada) frames were acquired at 30 Hz during a freehand scan, and exaggerated palpating movements were not necessary. The images are used only for qualitative assessment of the strain estimation algorithms. Results are shown for (1) olive/gelatin phantom mimicking a stiff inclusion in soft tissue, (2) tissue-mimicking phantom with two layers, (3) human male breast in vivo, (4) human calf muscle in vivo (see Figures 17-20 later).

For comparative purposes, we tested phase, correlation coefficient and adaptive strain estimators. The performance of phase and correlation coefficient estimators is compared for several variations: uncorrected strain estimation, log compression, limit log compression and AMC. Quantitative tests use simulation data, where the performance is measured by evaluating SNR₆ ; the strain standard deviation is calculated from the raw strain estimates, where no smoothing has been applied. For a qualitative assessment, we also present example images from in vitro and in vivo scans. Fair comparison is made possible by fixing the window parameters across all of the estimators in each test. It should be noted, that where there is a priori knowledge of a uniform strain field, the process of imaging strain by differencing closely spaced windows serves only to introduce noise; instead, windows separated by a large distance should be differenced in order to achieve an SNR _e that becomes arbitrarily high for large window spacing. Alternatively, in practical systems it is sensible to match larger numbers of closely spaced windows, and to combine their estimates by filtering methods such as least squares or wavelet decomposition. To varying degrees, these techniques reduce both noise and resolution, although the AM effect will remain important. To investigate the noise that is introduced by erroneous estimation location assumptions, and to evaluate the performance of the proposed AMC technique, in our quantitative tests we use the method of differencing windows at a fixed window- spacing, At = 2.7/1 (i.e. 0.45 μ s, 0.35 mm, 30 RF samples at 66.7 MHz). The window length, T , is varied between tests, with the chosen length stated in each case.

The experimental methods which follow provide a full description of each estimator, the properties of the simulation data, and the nature of the in vitro and in vivo ultrasound scans.

Efficient phase zero search

The efficient phase zero search (EPZS) adapts the concept of Pesavento et al. (ibid). To summarise, a 5-10 MHz filter is applied to the RF data (T₁ , r₂ ) before converting to analytic signal representations («, , a₂), which are modulated to the baseband (a_bl , a_b2) to enhance the accuracy of linear interpolation. a_b2 is estimated at subsample locations by baseband linear interpolation, to enable accurate subsample estimation of d (for a discussion of interpolation frequency responses, see Proakis and Manolakis (J. G. Proakis and D. G. Manolakis, Digital signal processing: principles, algorithms and Applications, Upper Saddle River, third edition, 1996). Phase-based methods begin with the displacement of the analysis window known already to within λ/2 ; this is achieved by initialising each window with the final displacement estimate from the preceding one; windows at the top of each A-line are initialised with d = 0 . Displacement estimates are differenced to produce strain estimates following Equation 1. The estimation location is usually assumed to be the window centre.

f_n = nAt + ~ (32)

The phase is preserved but the amplitude is partially suppressed when the signal is log compressed according to the following formula.

a_bM(_.t) = log(l + c \ a_b(t) \)e^Ml) (33)

c is the compression factor. The larger the value of c , the smaller the amount of amplitude information that is retained, since the size of variations in the log compressed amplitude becomes smaller compared to the mean value. We refer to this algorithm as EPZS_L1. As c — > oo all of the amplitude information is discarded, since log compressed amplitude variations become infinitely smaller than the mean. Limit log compression has a simpler form.

W⁾ = ^e""^tl° ⁽³4⁾

We refer to limit log compression as EPZS_L2. For phase-based methods, EPZS_L2 is the counterpart of one-bit compression or zero crossing techniques in correlation coefficient methods. We also present results for EPZS with AMC, referred to as EPZS_A. In addition to producing analytic signals, we detect the signal envelope, \ a \ , which is exploited as follows for AMC estimation of f „ (c.f. Equation 28). \ t

EPZS_L2 uses no amplitude information; so the AMC version of f „ following equation (35) can be identical to the window centre assumption. However, EPZSJLl still exhibits a degree of AM susceptibility, so results arc presented for an algorithm combining EPZS_L1 with AMC (operating on the log compressed signal envelope), referred to as EPZS_LA.

Correlation coefficient maximiser

The correlation coefficient maximiser (CCM) searches initially at integer sample locations for the maximum value of the cross correlation coefficient (see Equation 29). The estimate is refined by allowing subsample values of d and interpolating r₂ at subsample locations. Again, a complex baseband representation of r₂ allows highly accurate subsample interpolation, as with EPZS, but in CCM it is converted back to a subsample real signal for the correlation coefficient calculation. This requires the following calculation, where ω_m is the modulation frequency that was used earlier to shift the analytic signal down to the baseband.

r₂(t) = K{a_b2(t)e^jω»><} (36)

τ is again usually assumed to be the window centre (Equation 32). Log compression (CCM_L1) is tested as a means of reducing the error in f , using the following formula:

r_log (0 - log(l + c I r(0 \)sign(r{t)) (37)

To maximise algorithm performance, the full RF signal is used for subsample interpolation of r₂ , which is only log compressed at the moment of computing the correlation coefficient. In the limiting case as c → ∞ variations in the log compressed signal magnitude become infinitely smaller than the mean magnitude, so only the sign is important. A simpler expression may be used.

Subsample interpolation actually still employs the full RF signal, so zero crossings are identified with high accuracy. We call this variation CCM_L2, although it could be described as one-bit compression and is similar to a zero-crossing method. AMC is applied to CCM following Equation 31, which is referred to as CCM_A. AMC is also applied alongside non-limiting log compression in CCM_LA.

Adaptive strain estimator

Our adaptive strain estimators have two search dimensions — displacement and stretch — for each spatial dimension of strain estimation. The algorithm has the following stages: (1) each post-deformation window is shifted till the best match is located; (2) the shifted window is stretched to maximise a similarity measure; (3) displacement is re- estimated for the stretched window; and (4) the process repeats iteratively until convergence. Once arrays of displacement and stretch have been calculated, either the displacement estimates may be differenced (as in displacement-based methods) to re- estimate strain, or the stretch estimates may be displayed directly (which is the approach followed in this study). An estimator of this form was observed to produce significantly better strain images than those that are achieved by the basic displacement estimation approaches, with the greatest improvement for high strains. SAD was found to outperform the correlation coefficient, so this is the chosen signal similarity measure. The origin of this difference may lie in the fact that often p_rr — 1.0 at the correct stretch, in which case SAD is less prone to quantisation errors.

It has subsequently been noted that a minor modification to the adaptive stretching algorithm yields a further performance improvement. The modification concerns the way that displacement is estimated: our adaptive strain estimator (ASE) estimates the locations of the windows directly from the strain estimates, rather than searching over two dimensions. This has been found to yield higher SNR _e .

The initialisation of EPZS depends on the fact that the displacement at the top of each A-line is zero. Similarly, ASE searches only over strain (and not over displacement) in the top window of each A-line. This utilises the knowledge that a search over displacement could only degrade the accuracy of the estimate in the event that a nonzero displacement were found for the top of the window. The displacement at subsequent windows is estimated accurately by integrating the estimated strains, where it will be recognised that integration is a noise-suppressing operation. The offset of the first sample in a succeeding overlapping window is, of course, not equal to the displacement at the end of the first window. Rather, the relationship we assume is illustrated in Figure 7, where estimated strains are displayed as gradients on a plot of displacement against time. The window strain estimate multiplied by the window length, Ts,,-_{ι >} provides the best estimate for the displacement difference between the end and the start. The following window is therefore pinned at this end point, and stretched on either side to find the next estimate, S_n ■ This means that the offset displacement at the start of window n depends on: the offset of window n -l , the previous window stretch, and the candidate window stretch, S .

d_o.* = ^₁ + S-F-S(T - to) (39)

This leads immediately to a second result for increased resolution with overlapping windows. An estimate that resolves strain changes at the scale of the shift between windows (thereby matching the resolution of the displacement methods) is produced as follows.

£„_, = {$„_?- S_n(T -to)}/to (40)

This is a consequence of the geometry in Figure 7. Increased resolution comes at a cost of increased estimation noise. We present results using s rather than S , however, since the higher resolution of s makes it the appropriate comparison with the displacement- based methods.

SAD(n,s) = (41)

f is the sampling frequency. S_n minimises SAD(n,s) ,

S_n = arg min SAD(n,s) (42) It may be possible to adapt fast algorithms to this optimisation problem, but for now we use an exhaustive search.

Benefits of amplitude modulation correction

What follows is an analysis to show that, in general, discarding amplitude information is undesirable. We analyse a generic displacement estimator, motivated by the actual properties of phase-based methods. We assume that a window of arbitrary length produces an unbiased displacement estimate. The shortest possible window covers one RF sample, producing a displacement estimate, d' . The estimation variance, σ_ά ² , is inversely proportional to the local ultrasonic SNR; this assumption follows the Cramer- Rao (ibid) lower bound for displacement estimation variance.

We assume a simple model for RF signals during a strain imaging ultrasound scan. An underlying signal, r , is present in both the pre- and post-deformation signals, r_x and r₂ , but these are recorded in the presence of additive noise.

r, (O = KO + «,(') (44) r₂(t + d(t)) = r(t) + n₂(t + d(t)) (45)

H₁ and n₂ have zero mean, with power

. They are mutually uncorrelated, and both noise signals are uncorrelated with r . In general n_x and n₂ consist not only of electronic noise — other sources of uncorrelated signal components include morphological changes to the speckle pattern and non-axial scatterer motion. The SNR can be expressed in terms of these signal components.

The constant of proportionality in Equation 44, C₁ , must be a large number, since the short windows produce inaccurate estimates. However, the generic estimator actually uses longer windows, yielding a weighted average of the single-sample estimates.

d_n is the final displacement estimate at window n , and W(t) is the weighting for estimate d'(t) . If errors in the single-sample estimates are mutually uncorrelated, then the variance of the overall estimate is as follows.

This can be minimised by choosing W as follows, where C₂ is an arbitrary constant.

C₂ C₂SNR(J.)

W(t) = ^- = ^^> (49) σΛ² t) C

The implications of this result are not immediately obvious, since SNR (t) is an unknown quantity. However, the expected error is minimised by choosing W according to the expected value of the local SNR, given the information that is available. We require the statistical expectation of the RHS in Equation 47.

= E[r{tf ]xE[2{n_l{tf + n₁{t + d(t))²r] (51)

The expected noise term is assumed constant ( C₃ ). More sophisticated noise estimates are possible if assumptions can be made about the statistical properties of the noise source, but we restrict ourselves to the most general approach (note, E[V¹ j ≠ E[^]^"1 , so

C₃ ≠ σ;² ).

E[SNR(t)] = C₃E[r(tf ] (52)

Since the noise is uncorrelated and the displacement estimate is assumed to be unbiased, the expectation of the local cross power of the recorded signals is equal to the expected signal power.

E[r_l(t)r₂(t + l)] = E[r_l(t)r₂(t + d(t))] (53)

= E [(KO + Ti₁ (O)(KO + H₂ (t + d(t)))] (54)

= E[r{tf ] + E[KO^₁(O] + E[r(t)n₂(t + d(t))] + E[n_{{t)n₂{t + J(O)] (55)

= E[K0²] (56)

The cross power can therefore be taken as an estimate of the signal power. By combining the results of Equations 50, 53 and 57, it emerges that the optimal weighting for each single-sample displacement estimate can be evaluated. In the following expression C₄ is an arbitrary constant.

W{t) = C_Λ(t)r₂(t + d_β) (57)

Weighting by this formula minimises the expected value of σ²- .

Experimental results

Quantitative results indicate the advantages and disadvantages of each technique. Important trends are illustrated by graphs. Where there is space for error bars these extend to one standard deviation either side of the mean. We also present strain images for qualitative assessment.

Window lenfith

Results for EPZS, EPZS_A, CCM, CCM_A and ASE with window lengths, T , in the range 2.8-27.1 λ indicate a suitable choice of T for the later tests. They also serve as a first opportunity for assessing the AMC technique. Figures 8 and 9a show performance against window length at 0.5% strain, while Figure 9b shows the effect of window length on EPZS_A and ASE at a higher strain. 13.52 is employed for all other results herein.

More particularly Figure 8 shows SNR _e against window length for EPZS and EPZS__A, with both 71 dB and 20 dB data at 0.5% strain. Uncorrected EPZS with 71 dB data reaches a plateau at T = 10/1 , which the 20 dB results converge towards for long windows. When AMC is applied there is no such plateau and much higher SNR _e is achieved — SNR _e is initially a linear function of window length, and it continues to increase for long windows, although the gradient becomes less steep.

Figure 9a shows SNR _e against window length for CCM and CCM_A, with both 71 dB and 20 dB data at 0.5% strain. Uncorrected CCM is almost identical to EPZS. However, AMC is obviously less accurate for CCM, since the improvement with CCM_A is much smaller and the results are erratic for long windows. Figure 9b shows SNR _e against window length for EPZS_A and ASE, with 20 dB data at 4% strain. ASE performs less well with short windows, but it reaches a high and fairly constant level of performance for T > 10/t . EPZS_A, by contrast, performs well with short windows and has a higher peak SNR _e . However, windows with T > lθλ have a differential displacement of

>0.4/t between the ends, so in this range EPZS_A suffers substantially increased decorrelation and estimation noise.

To illustrate the practical meaning of SNR₂ , Figure 10 shows strain images at 0.5% compression. More particularly, Figure 10 shows strain images for a 0.5% compression with 20 dB data using T = \3.5λ : (a) EPZS; (b) CCM; (c) EPZS_A; (d) CCM_A. The performance of EPZS and CCM is similar, though EPZS A performs considerably better than CCM_A.

The characteristics of the images can be compared with the corresponding SNR₆ results from the graphs. The images have a linear scale with 0 (black) representing zero strain, 127.5 (mid-grey) is the simulated strain of 0.5% and 255 (white) represents 1%. Saturation occurs at 0 and 255, and no smoothing has been applied, so each section between successive estimation locations has constant brightness. An ideal estimator would yield a uniform greyscale level, but this is unachievable in practice.

Compression factor

A justification is presented for the choice of log compression factor in the algorithms EPZS_L1, EPZS_LA, CCM_L1 and CCMJLA. The effect of log compression varies to a large degree depending on the strain level, so Figures 11-13 show results at strains representing the smallest, largest, and mid-range in the simulation data. It is evident that log compression is not always desirable, but the choice of c reflects a value that is likely to boost SNR _e in high strain regions, whilst avoiding extreme degradation of low strain estimates, c = 10³ is employed for all of the remaining results.

More particularly, Figure 11 shows SNR₆ results for EPZS-Ll, EPZS_LA, CCM_L1 and CCM_LA with 20 dB data at 0.01% strain as a function of c , the compression factor. At low strains, the main effect of log compression is increased noise. This effect is more pronounced with CCMJLl. AMC has almost no effect in these images. Figure 12 shows SNR _e results for EPZS-Ll, EPZSJLA, CCM_L1 and CCM_LA with

20 dB data at 0.5% strain as a function of c , the compression factor. At this strain, log compression significantly improves the performance of EPZS_L1. CCM_L1 is also improved by slight log compression. Better performance is produced by AMC, although this is degraded by log compression, so as c — > ∞ EPZS_LA and CCM_LA converge with the curves where AMC has not been applied. Figure 13 shows SNR _e results for

EPZS-Ll, EPZS_LA, CCM_L1 and CCM_LA with 20 dB data at 4% strain as a function of c , the compression factor. At this strain all of the algorithms can be improved by applying an appropriate level of log compression. The greatest improvement is exhibited by EPZS-Ll, while the ACM algorithms are still degraded by high compression factors, and they eventually converge with the curves where AMC has not been applied. Strain dependence

With parameters T and c selected as per the preceding sections, Figures 14-16 compare the performance of EPZS, EPZS-Ll, EPZS_L2, EPZS_LA, EPZS_A, CCM, CCMJLl, CCM JL2, CCM_LA, CCM_A and ASE across a range of strains. More particularly, Figure 14 shows SNR _e -strain characteristics for the EPZS family of algorithms with 20 dB data. EPZS_A has the best performance across a wide range of strains, although the SNR₆ is lower at high strains and at 4% the best performance is from EPZS_LA. Figure 15 shows SNR₆ -strain characteristics for the CCM family of algorithms with 20 dB data. At all strains CCM_A significantly outperfoπns the other algorithms. In the absence of AMC, log compression boosts CCM performance at strains above 1.5%, though the best log compression performance comes from the combination algorithm, CCMJLA. Figure 16 shows SNR _e -strain characteristics for

EPZS_A, CCM_A and ASE with 20 dB data. These are the best algorithms from each of the three families. EPZS_A performs best across most strains, though ASE does slightly better at 4%, where the other algorithms have lower SNR _e owing to significant decorrelation.

In vitro and in vivo results

Finally, images from real ultrasound scans are presented in Figures 17 to 20. For the sake of conciseness, we restrict ourselves to EPZS, EPZS-Ll, EPZS JL2 and EPZS_A, allowing a qualitative assessment of log compression and AMC when applied to real data. The images in Figures 17-20 have been smoothed slightly by estimating strain with a 1 mm least squares filter along the axial direction; no other filtering has been applied and the values of parameters T and c are unchanged.

Interpretation of results

Window length results in Figure 8 show that AMC is extremely effective when applied to EPZS, which validates the above analysis. Notice that while increasing the window length is known to reduce σ², , nevertheless the unconⁿected algorithm quickly reaches a plateau: this is because the primary source of noise is the AM effect when long windows are used. Meanwhile, when AMC is applied the remaining noise is mainly due to σ²- , so higher SNR _e is achieved with the 71 dB data. It is encouraging, however, that the curve for 20 dB data has the same form as for 71 dB data. This shows that although AMC was derived considering noiseless data, the technique has a similar effect in the presence of noise.

Note from Figure 9a that the performance of uncorrected CCM is almost identical to EPZS. However, AMC for CCM is less successful, which probably reflects the lack of a formal derivation, rather than implying that it is not possible to correct the AM effect in this case. The formula in Equation 31 was based on intuition based on our previous insights since the derivation of a superior CCM_A algorithm is a considerably more challenging mathematical problem than EPZS_A.

Figure 9b confirms that ASE offers an alternative route to high-performance strain estimation, hi particular, it is possible to achieve good performance using arbitrarily long windows. This means that locations of extremely high strain will not be subject to reduced SNR _e when the window length has been chosen for optimal performance at a range of lower expected strains. It is also interesting to note that EPZS_A actually outperforms ASE for short window lengths, and EPZS_A has the higher peak performance. EPZS_A performs less well with longer windows, where high strains cause significant decorrelation. The window length chosen for subsequent tests (T = 13.5/1 ) was determined by two factors: (1) long windows eventually reduce resolution in practical strain imaging; and (2) 13.5/1 is a sensible balance for near- optimal performance across all algorithms at all strains in the range 0.01-4%.

The log compression results in Figures 11-13 demonstrate the behaviour that we predicted in the discussion of the benefits of AM above. σ_δ ² dominates at the low strain in Figure 11, so preferably noise suppression uses all of the amplitude data to maximise the accuracy of the displacement estimates. Therefore, log compression degrades performance. EPZS and EPZS_A have identical

, while σ| is less important, so

AMC is of less benefit. The same observation applies to CCM and CCM_A. However, EPZS and EPZS_A are degraded less severely by log compression, since the retention of phase information makes these algorithms more robust. CCM only uses the real signal, so σ_& ² increases rapidly with log compression as information is discarded.

However, 0.5% strain in Figure 12 is already sufficiently high for the noise contribution of σ_p ² to become important. Log compression yields a significant improvement in

EPZS_L1, and slight log compression also improves CCM_L1. Better performance is achieved by the AM corrected algorithms, although these are still degraded by log compression. EPZS_A and CCM_A eventually converge with the uncorrected curves as c -> co . Log compression is most beneficial at the higher strain in Figure 13. Estimation noise here comes mostly from σ? , so EPZSJLl performs much better when a high level of log compression is applied. CCM_L1 is also improved by high log compression, although it peaks at a relatively low value of c . The AMC algorithms are also improved by slight log compression, indicating that the AMC formulae are less accurate at high strain, so a combination of AMC and log compression yields the lowest location variance. However, the AMC algorithms have considerably higher peaks than the uncorrected algorithms, so performance convergence as c -> ∞ represents a significant reduction in SNR _e . The choice of c = 10³ for subsequent tests reflects a balance between the EPZS_L1 optima at 0.5% and 4% strain.

The SNR _e -strain characteristics in Figures 14-16 further demonstrate the advantage of applying AMC. It yields the best performance in both EPZS and CCM families of estimators. The uncorrected EPZS and CCM curves again reach a plateau in the region where σ², dominates, as predicted by Equation 7. It is also interesting to note that the AMC curves peak at lower strains than the other algorithms, which follows from the combined effects of AMC becoming less accurate at high strains and σ? becoming more important as the level of signal decorrelation increases. In the case of EPZS_A, AMC is precisely accurate for small strains, but it diverges from the correct estimation location at higher strains where errors in the assumptions of the derivation become increasingly significant. The hybrid algorithm, EPZS LA is the best at 4% strain, so the combination of AMC with moderate log compression may be the best noise suppression strategy at high strains. Figure 16 compares the best estimators from each family of algorithms. EPZS_A has the best performance at most strains by a large margin. At low strains the worst algorithm is ASE. This may indicate that the signal stretching technique is inherently more noisy, although at higher strains its advantages are the absence of the AM effect and lower signal decorrelation. Therefore, ASE outperforms CCM_A for strain >2%, at 4% it also outperforms EPZS_A, and the gradient of the curve is still positive, so ASE may offer further performance benefits at yet higher strains. However, it is likely that the main advantage of ASE is the relative independence of performance and window length. On the other hand, we have already seen in Figure 9a that EPZS_A can (depending on window length) outperform ASE by a large margin.

Images from real ultrasound scans in Figures 17-20 provide further evidence of the comparative properties of these algorithms. In general, the EPZS_A images are the least noisy, while EPZS_L1 and EPZS_L2 are more or less noisy than EPZS depending on the local strain (c.f. Figures 11-12). These images also demonstrate the importance of AM correction when AM artefacts correlate with features in the B-scans. It is evident in Figure 17 that the AM effect has distorted the shapes of features in the EPZS image, particularly in the attenuation shadow below the olive. Figure 18 shows a more extreme example. The specular reflection is of unknown origin — possibly a crack has developed in the gelatin matrix. It causes severe distortion of EPZS, where the dark patch in Figure 18b looks like a low strain planar inclusion. However, EPZSJL2 is provably unaffected by amplitude variations, so real tissue features must also appear in Figure 18d. The dark patch is absent, proving that it is actually an artefact. A mild artefact is also observed with EPZS_A in Figure 18e, where the local sparseness of estimation locations around the reflection causes a textural change in its vicinity.

The in vivo images in Figures 19 and 20 demonstrate that AM artefacts often occur in scans of real human tissue — isoechoicity is rarely a feature of salient scan planes. The male breast in Figure 19 has an appreciably different strain image with EPZS compared to the other algorithms. A bright band at the top of the lower section reappears as a zero- strain band in Figure 19b, but this is an artefact, absent from Figures 19c-e. Many similar artefacts are present in the calf scan of Figure 20. This is extremely anisoechoic, and comparison between Figure 20b and Figures 20c-e shows that all of the main features in the EPZS image are artefacts.

To summarise, the AM effect has been theoretically introduced and empirically investigated and a new technique we call AMC has been derived for the enhancement of ultrasonic displacement and strain estimates. Simulation, in vitro and in vivo results show a substantial reduction in the level of estimation noise. However, it is always possible to reduce noise by applying filters, thereby sacrificing spatial resolution in order to boost SNR _e . It is likely in practice, therefore, that the main impact of AMC will be an improvement in spatial resolution, and AMC can be used to enhance strain imaging in 2D or 3D if required. The ultimate limiting factor in ultrasonic displacement and strain estimation will be the limited bandwidth of RF ultrasound signals, i.e. the point spread function is not an impulse. This means that even if signal displacements were tracked perfectly, there would be a residual error between those displacements and the actual tissue motion. Developments in ultrasound deconvolution for enhanced ultrasonic resolution may help on this point, hi principle, AMC has other more general applications for enhanced tracking of small motions.

When AMC is applied with regularly spaced windows of a fixed length this leads to variable spacing of the estimation locations. It may therefore be helpful to automatically vary the length and spacing of the windows to maintain spatial resolution with AMC, or to achieve a balance between spatial resolution and estimation noise according to an appropriate cost function.

We have used an assumption of locally constant strain and estimation noise will increase when the second derivative of displacement is non-zero within any particular window. These first order corrections are already very useful, but potentially superior AMC formulae could be achieved by exploiting correlations between the errors in overlapping windows. Log compression in both its moderate and limiting forms has been demonstrated to be particularly useful with phase (EPZS_L1 and EPZS_L2) and the retention of phase information, regardless of how far the amplitude is compressed, makes phase-based methods more robust. In some tests EPZS_L2 has been one of the most successful algorithms. This may be because EPZS_L2 exhibits a high level of displacement estimation error (σ|) and EPZS_L2 is unaffected by location errors (σj, ).

At high strains σjj, is often the primary source of error, so EPZS_L2 can outperform some of the other estimators and there are computational advantages if all of the amplitude information can be discarded. However, it is inferior to the EPZS_A algorithm incorporating AMC, which is marginally more computationally expensive, but still suitable for real-time strain imaging. EPZS_A outperformed all of the other algorithms tested in this study throughout the typical range of strains encountered in practical ultrasonic strain imaging systems.

No doubt many other effective alternatives will occur to the skilled person. It will be understood that the invention is not limited to the described embodiments and encompasses modifications apparent to those skilled in the art lying within the spirit and scope of the claims appended hereto.

Claims

CLAIMS:

1. A method of processing at least one-dimensional image data captured by an imaging technique to determine deformation data defining an at least one-dimensional deformation in an imaged object, said deformation data comprising at least one data pair, said data pair comprising displacement data and corresponding displacement location data, the method comprising: inputting first and second sets of said image data, a said set of image data comprising imaging signal data including at least one of signal magnitude data and signal phase data for an imaged region of said object; positioning a first window at a first position on said first set of image data and determining a corresponding position for said first window in said second set of image data; determining first displacement data defining a first displacement estimate from said positions of said first window in said first and second image data; and determining an estimated location of said first displacement estimate responsive to said imaging signal data within said first window for at least one of said first and second sets of image data; whereby said at least one data pair comprises said first displacement data and said estimated location of said first displacement estimate.

2. A method as claimed in claim 1 wherein said deformation data comprises a set of said data pairs, the method further comprising: positioning a second window at a second position on said first set of image data and determining a corresponding position for said second window in said second set of image data; determining second displacement data defining a second displacement estimate from said positions of said second window in said first and second image data; and determining an estimated location of said second displacement estimate responsive to said imaging signal data within said second window for at least one of said first and second sets of image data; whereby a second data pair of said set of data pairs comprises said second displacement data and said estimated location of said second displacement estimate.

3. A method as claimed in claim 1 or 2 wherein said imaging signal data includes said signal magnitude data, and wherein said determining of an estimated displacement estimate location for displacement data determined from a said window is responsive to said signal magnitude data within the window.

4. A method as claimed in claim 1, 2 or 3 wherein said imaging signal data includes said signal magnitude data, and wherein said determining of an estimated displacement estimate location for displacement data determined from a said window is responsive to said signal magnitude data, for the window position, for both said first and second sets of image data.

5. A method as claimed in claim 4 wherein said determining of an estimated displacement estimate location for displacement data determined from a said window comprises determining a centroid using said signal magnitude data for the window position for one or both of said first and second sets of image data.

6. A method as claimed in any preceding claim wherein said imaging signal data includes said signal phase data, and wherein said determining of an estimated displacement estimate location for displacement data determined from a said window is responsive to said signal phase data for one or both of said first and second sets of image data within the window.

7. A method as claimed in any preceding claim wherein said determining of an estimated displacement estimate location for displacement data determined from a said window uses a weighting function determined responsive to a procedure for said window positioning.

8. A method as claimed in any preceding claim wherein said each set of image data comprises a plurality of digitised image data samples, and wherein said window position determining includes interpolating between positions corresponding to said data samples.

9. A method as claimed in any one of claims 1 to 8 wherein said corresponding window position determining comprises cross-correlating image data from said first and second sets of image data.

10. A method as claimed in claim 9 wherein said corresponding window position determining further comprises determining a zero-crossing of a phase of said cross- correlation.

11. A method as claimed in claim 9 wherein said corresponding window position determining further comprises determining a maximum correlation coefficient of said cross-correlation.

12. A method as claimed in any preceding claim further comprising positioning a further window at a further position on said first set of image data and determining a corresponding position for said further window in said second set of image data; and determining displacement data defining a further displacement estimate for said further window and an estimated location of said further displacement estimate to provide a further said data pair for said deformation data.

13. A method as claimed in any preceding claim wherein said image data comprises two-dimensional data, and wherein said deformation data comprises a two-dimensional deformation data.

14. A method as claimed in any preceding claim wherein said deformation data defines a strain field within said object.

15. A method as claimed in any preceding claim further comprising determining one or more strain values of said object from said deformation data.

16. A method as claimed in any preceding claim further comprising inputting stress data defining a stress applied to said object, and determining elasticity data for said object from said stress data and said deformation data.

17. A method as claimed in any preceding claim further comprising displaying an image of said deformation data.

18. A method as claimed in any one of claims 1 to 17 wherein said imaging technique comprises a pulse-echo imaging technique, and wherein said imaging signal data comprises pulse-echo signal data.

19. A method as claimed in claim 18 wherein said pulse-echo technique comprises ultrasound imaging.

20. A method as claimed in claim 18 wherein said pulse-echo technique comprises magnetic resonance imaging.

21. A method as claimed in any preceding claim wherein said object comprises biological tissue.

22. A carrier carrying processor control code to, when running, implement the method of any preceding claim.

23. Apparatus for processing at least one-dimensional image data captured by an imaging technique to determine deformation data defining an at least one-dimensional deformation in an imaged object, said deformation data comprising at least one data pair, said data pair comprising displacement data and corresponding displacement location data, the apparatus comprising: an input for first and second sets of said image data, a said set of image data comprising imaging signal data including at least one of signal magnitude data and signal phase data for an imaged region of said object; a system for positioning a first window at a first position on said first set of image data and determining a corresponding position for said first window in said second set of image data; a system for determining first displacement data defining a first displacement estimate from said positions of said first window in said first and second image data; and a system for determining an estimated location of said first displacement estimate responsive to said imaging signal data within said first window for at least one of said first and second sets of image data; and wherein said at least one data pair comprises said first displacement data and said estimated location of said first displacement estimate.

24. Apparatus as claimed in claim 23 wherein said deformation data comprises a set of said data pairs, the apparatus further comprising: a system for positioning a second window at a second position on said first set of image data and determining a corresponding position for said second window in said second set of image data; a system for determining second displacement data defining a second displacement estimate from said positions of said second window in said first and second image data; and a system for determining an estimated location of said second displacement estimate responsive to said imaging signal data within said second window for at least one of said first and second sets of image data; and wherein a second data pair of said set of data pairs comprises said second displacement data and said estimated location of said second displacement estimate.

25. Apparatus as claimed in claim 23 or 24 wherein said imaging technique comprises a pulse-echo imaging technique, and wherein said imaging signal data comprises pulse-echo signal data.