Det medisinske fakultet

Mathematics of strain and strain rate.

Edited by Asbjørn Støylen, dr. med.


This section updated: November 2010

This section is dedicated to a more in depth treatment of the mathematical basis of strain and strain rate imaging, for those with special interests. Mathematicians and engineers, however, will find this embarrassingly simple, and should rather refer to the much more in depth analyses in (Andreas Heimdal. Doppler based ultrasound imaging methods for non-invasive assessment of tissue viability, NTNU 1999.). On the other hand, the present section does not use advanced mathematics, out over what is obtained by a general education.

Section index


Other sections:

Basic principles of ultrasound and scanner technology.

Is deformation imaging useful?

Measurements of strain and strain rate by ultrasound

Back to main website index

The Doppler effect.

The Doppler effect is the effect of the velocity of the observer (A) or the wave source (B) on the perceived wavelength, as originally derived by CA Doppler (119). The basic fact is that the velocity of a wave, c, is constant in a given medium, and equal to the number of oscillations per second, times the wavelength (length of one oscillation):
 
The time of one oscillation, i.e. the time it takes for the wave to move one wavelength is then:


A
 B

The moving source is illustrated in an animation in: http://folk.ntnu.no/stoylen/strainrate/Ultrasound/index.html#Doppler

Thanks to Hon Chen Eng of University of Toledo who pointed out an inconsistency in the original illustration
and showed a better way of illustrating the Doppler effect in this image.
An observer (blue) moving towards a stationary wavesource with the velocity:

will meet the wave as the wave have moved a distance , which is the perceived wavelength. The observer has moved the distance:

The motion of the wave and the motion of the observer happen during the same time interval:

The motion of the observer thus shortens the original wavelength  by , i.e:

Thus:



The change in frequency, the Doppler shift is:
If the source moves toward a stationary observer with the velocity:             

In the time the wave moves one wavelength, the source moves the distance:                                                                            

The motion of the wave and the motion of the source happen during the same time interval:

The distance from the next wave emitted from the new position of the source (small dotted red circle)
to the observer (blue) is shortened by  in the direction of the motion, so the new wavelength representing the distance between the first and second waves is

Thus:

The change in frequency, the Doppler shift is:




If  v << c, then:


and:


For reflected ultrasound, the effect is twice as great. A reflector moving towards the source will shorten the incoming wavelength in the same way as an observer moving towards the source, and the reflected ultrasound wavelength will be further shortened in the same way as a moving source following the reflected ultrasound, Thus, the effect is:

(The approximation in case B is small, the velocity of ultrasound in tissue is 1540 m/s, while the velocity of blood is between 0,2 and 6 m/s, and tissue between 0.05 and 0.2 m/s, giving a v/c of maximum 0.004, i.e the approximation in B is maximum 0.4% and in reflected ultrasound 0.2%). In all cases, if the velocity vector has an angle with the direction of the observation, the effective velocity is v  as illustrated below, in the case of reflected ultrasound, the angle is the angle between the emitted ultrasound beam and the  velocity vector. Thus the full Doppler equation for reflected ultrasound is:

In all cases it is evident that velocities at right angle to the ultrasound beam will result in no Doppler shift, and if the reflector moves away from the ultrasound source, there is a negative Doppler shift.


Phase analysis

In colour Doppler, the analysis is done in terms of the phase shift: A wave can be described as a sine wave, and thus, any point on the wave can be described by which phase of the wave the point is, as illustrated below:

Phase analysis.  If the waveform is treated as a sine curve,  every point on the curve corresponds to an angle, and the phase of the point in the curve can be described by this angle; the phase angle .  From the diagram, it's also evident that a full wavelength, , is equivalent to 2, and for every point the corresponding fraction of a wacelength is equivalent to an angle which is the fraction of 2. However, from the diagram at the top, it is evident that by sampling the waveform only once, the phase is ambiguous, it is not possible to separate the phase of point a from point b.  The two points are separated by a quarter of a wavelentgh, or 90° (). In order to determine the phase of the points unambiguously, the pulse has to be sampled at to points separated by less than a quarter wavelength. Then it can be seen that point a is in increasing phase from a1 to a2,  corresponding to a phase angle of  0 - /2 while b is in a decreasing phase corresponding to an angle of /2 - .

Shooting at least two (or more) pulses in rapid sequence, (NOT to be confused by sampling one pulse at two timepoints as illustrated above) results in the possibility to analyse the Doppler shift in terms of the phase shift between the pulses. The phase shift analysis is based on the principle that when pulse 2 hits a moving scatterer, the scatterer will have moved a little away from, or towards the probe, and the return pulse 2 will then be in  a different phase from pulse 1. The distance the scatterer has moved, is of course a function of the velocity of the scatterer and the time between pulses given by the pulse repetition frequency, thus:
d = v * t = v * 1/PRF.






Two pulses sent toward a scatterer with a time delay  t2 - t1 = 1/PRF. Given that the scatterer has a velocity, it will have moved a distance, d, that is a function of the velocity and the time (d = v x t).  Thus, pulse 2 travels a longer (or shorter) distance equal to d with the speed of sound, c, before it is reflected.  During the time pulse2 has travelled the distance d to the new position of the scatterer and back to the point of the reflection of pulse 1, i.e. a distance 2d  pulse 1 has travelled the same distance away from the reflection point. (The scatterer will have travelled further, but this is not relevant).  Thus the diasplacement of the waveform of pulse 2 relative to pulse 1, is 2d. This corresponds to a phase shift from pulse 1 to pulse 2 of  , and it can be shown that = 4  f0 / v PRF By sampling the two pulses simultaneously at two timepoints, as shown in the previous illustration, the phase of each pulse can be determined. The analysis  f hte relative positions of all four points is done by autocorrelation, a quick (and dirty?) method that allows online computation.






As descibed above, a pulse has a certein bandwidth, describing the frequency content of the pulse. In spectral analysis, this will give a spectrum of a certain width, corresponding to the velocity distribution of flow velocities. In phase analysis, this will correspont to a certain distribution of phase angles as illustrated. Autocorrelation, however, will only result in the average phase angle.
 In the case of stationary noise (clutter) as f.i. reverberations, the autocorrelation will result in an average phase angle that is in between the signal and the noise. The clutter noise will have to be removed by a low velocity filter in order to avoid severe underestimation of flow velocities.




In order to avoid aliasing, the distance the scatterer moves should be less than half a wavelength, i.e. the Nykvist limit. This is thus given by the pulse repetition frequency, exactly as in pulsed Doppler. In this case, the PRF is 1/t, the time delay between the two pulses in one package. In order to obtain an acceptable Nykvist limit, the PRF need to be at least 1 KHz which will result in a Nykvist limit of 19.8 at a depth of 15 cm. 


Thus, as a wavelength equals 2, the displacement (2d) of pulse 2 in relation to pulse 1 relative to equals the phase shift relative to 2:    2d/= /2. As d = v*t, t = 1/PRF and = c/f, it follows that:

= 4  *f0 / v*PRF

Thus, the phase analysis is based on the phase shoft from one pulse to next, instead of the apparent shortening of the wavelength by the motion of the source or observer (or reflector, which amounts to the combination of both. However, the phase shift from one pulse to the next is completely equivalent to the change in wavelength in one pulse, so the two are equivalent for all practical purposes.

The motion of the scatterer, d, is the same as the motion in the derivation of the Doppler equation, and the  displacement between the pulses, 2d is the same as the combined Doppler equation for reflected ultrasound: the displacement of sorce/observer as a fraction of the wavelength.
 




Lagrangian vs. natural and Eulerian  strain:

The Lagrangian definition of strain is:


This gives an intuitive approach, as positive strain is stretching, negative strain is shortening. But this is only relative to the original length. When an object is stretched from L0, and then is compressed, the strain is still positive till it is compressed beyond Lo, and vice versa, if it is initially compressed and then stretched, Lagrangian strain remains negative till it is stretched past L0. It may be more intuitive to define the instantaneous strain. As:

then:

Lagrangian strain is unknown if the original length is unknown, even if the instantaneous length and strain rate is known. Thus, a more natural approach may be to define instantaneous strain in relation  to the instantaneous length (Eulerian instead of Lagrangian coordinates):

The natural (Eulerian) strain is then:

The incremental increase in length is:

Summing all increments from t = 0 to t gives L(t+dt+dt) - L(t+dt) + L(t+dt) -L(t), and as L(t) = L0 at t=0 and L(t) = L at t, this gives the Eulerian strain:

Thus:

But as:

Then:

Thus Eulerian and Lagrangian strain are related:

and


The strain rate is strain per time unit, i.e. the temporal derivative of strain:



Velocity gradient and strain rate:

The velocity gradient was originally defined as the transmural velocity difference across the wall, divided by the instantaneous wall thickness W:

But this can be generalised to a general definition in all directions; the velocity difference between two material points, divided by the instantaneous difference between them:

The distance L changes with time, if v1 and v2 are different. The unit of the velocity gradient is cm/s/cm, which is equal to s-1. The gradient was found as the slope of the linear regression of the tissue velocities. (The linear regression assumes that the velocity distribution is homogeneous.) under the assumption that the velocity gradient is constant over the length L (spatially homogeneous).  By this, the longitudinal velocity gradient is:

But as:

then:
and:

In other words, the velocity gradient equals the rate of Eulerian strain and the Eulerian strain equals the time integral of the velocity gradient.

Relation to tissue Doppler:

In colour tissue Doppler, the velocities of all pixels are sampled simultaneously. The strain rate is sampled as the velocity gradient with a fixed offset distance x. As there is no tracking of the endpoints of the initial length (object), neither the velocities nor the distance relates directly to the strain of a defined object. We do not measure the instantaneous length change, nor the velocities at the endpoints. The strain rate estimator is thus:


This is not identical to the velocity gradient defined above, as L changes, but Dx is constant, except at the time when Dx equals L, then v(x) = v2 and v(x+Dx) = v1. However, Usually, however, L will differ from Dx, for most frames and objects, and the velocities will hence differ too. Under the assumption that the strain is equally distributed over the length of the object (spatially constant), SR will still be equal to the velocity gradient, i e the values of the two ratios will be the same. Strain being spatially constant means that the velocity increases linearly along the length as shown in the diagram:

For any L that is different from Dx, v2 – v1 will be greater or smaller than v(x) – v(x + Dx) by the same ratio. In the figure, this is evident, as the slope of the curve is the same wherever it is measured. As v1 and v2 are the velocities of the end points of L, the ratios SR and VG will be the same:

and the strain rate by tissue Doppler (SR) equals Eulerian strain rate.

Practical consequences of the difference between Eulerian and lagrangian strain

The difference between Langrangian and Eulerian strain rate has practical consequenses. When you measure the Lagrangian strain rate, it is the velocity difference relative to a fixed length, the initial length. In systole (where the length is decreasing), you measure the Eulerian strain rate, it is the velocity difference relative to a decreasing length, the instantaneous length. Thus, the Eulerian strain rate will be higher in magnitude than the Lagrangian, and the difference increases with increasing magnitude. In diastole, the differences will decrease, as long as the Lagrangian strain is calculated from the initial end diastolic (start systolic length). If, however, the end systolic length is taken as initial length, the instantaneous length is increasing, and the  Euleriean strain rate will be lower than  the Lagrangian.

The practical consequenses are that not only will peak systolic strain rate be higher in magnitude by Eulerian strain rate, but the timing will be later as shown below.
Left: two systolic velocity curves with velocities somewhat arbitrarily chosen, but within normal range, and with a fairly normal shape.  The values on the cyan curve is two thirds of the red one, and both have peak velocity at 0.05 sec.  Left: strain rate curves calculated from the velocity curves on the right side.  Assuming an initial distance between the points of the two velocity curves of 2.5 cm, again a realistic distance for the velocity difference.  Lagrangian strain (yellow) is derived by dividing the velocity difference at each time point by the initial distance.  This results in a simple inversion of the velocity difference, in a different scale, and with minimum value (maximum absolute) at 0.05 sec. Eulerian strain (green curve) is obtained by dividing with a strain length that decreases by the time interval times the velocity difference. In this case the differences between the curves is evident, and the minimum value (maximum absolute) of Eulerian strain rate is reached at 0.25 seconds.




This hase even more consequenses whan strain is derived from integration of strain rate, because the integration process will add the differences between the two measurements, as each little difference is added to the integral, much like interests on a capital.










Strain measured by integration of the strain rate curves in the figure above. It is evident
that the total systolic strain is about 6% higher (absolute values) by Eulerian strain rate. 		Strain as a function of integrated strain rate. Integration of velocity gradient leads to Eulerian strain (dotted line: identity). Lagrangian strain will give values that lies above Eulerian strain, i.e. negative values are lower in absolute values, while positive values are higher as one integrates from zero. (After Hans Torp by permission).
 Strain curves by the two methods: Grey as they appear by direct integration of the velocity gradient (strain rate), i.e. Eulerian or "natural" strain.  Black the Lagrangian strain by appying the correction formula below. (After Hans Torp by permission).
Thanks to Eirik Nestaas who pointed out that the colours of Lagrangian and Eulerian strain had been switched in the figure, this is now corrected.

Strain by integrated strain rate:

It has been shown that the natural strain rate can be estimated by tissue Doppler:

but then:

In commercial tissue Doppler, this correction is already applied, i.e. The scanner measures natural strain rate and integrates to Lagrangian strain. If strain is measured directly by other methods, for instance speckle tracking in ultrasound grey scale images or by MR as Lagrangian strain, in deriving strain rate, the inverse correction has to be applied to compare with standard strain rate by tissue Doppler (SR):



or simply:
 

As can be seen from the equation, the fomula converts the strain rate at a specific time from Lagrangian to Eulerian strain rate, but as shown above, the peak strain rate is not simultaneous by the two calculations. Thus, peak Lagrangian strain rate (which is lower than Eulerian), will be converted to Eulerian strain at the time of peak Lagrangian strain, not the peak Eulerian strain itself, so the conversion is practically useless. This is in contrast to the conversion between strains, as bot Eulerian and Lagrangian strain are cumulated maximal strrain, which is at end systole and simultaneous.

Strain by speckle tracking:




Speckle tracking can be done by a two-dimensional search. Defining a kernel region in the myocardium will define the speckle pattern within. The initial frame is shown in red. Within a defined search area (marked in blue), the new position of  this kernel can be recognised by finding the same speckle pattern within a like-sized frame in a new position. This indicates that each speckle has moved the same distance in the same direction (thin blue arrows), and the movement of the whole kernel then will have been the same (thick blue arrow). The size of the kernel defines the spatial resolution, and the size of the search region is defined by the maximal expected displacement from frame to frame. Higher frame rate will mean a smaller search region, if the velocity is the same. In practice, the speckle pattern does not repeat perfectly. However, every kernel has a unique speckle pattern due to the random nature of speckles. Finding the new position of the kernel can then be reduced to finding the like-sized area with the smallest difference or error in total pixel intensity with a trial matching kernel sized region within the search area. This is called the sum of absolute differences (SAD).

Where K is the original kernel area and Kt is a like sized area in the new location. The new kernel position is the area with the smallest SAD within the search region. This has been shown to track well-developed speckle patterns as accurately as normalized cross correlation (26).

In practice and in two dimensions, the algortihm works as (27):
Cross correlation can be used to weight the movement of the original kernel region to the kernel region with the lowest value for SAD method (28):



Strain in three dimensions

See also the section on strain in strain rate imaging.
Strain in real objects occur in three dimensions. This complicates the matter, as the strain tensor then has more components, the number of components increase by the square of the number of dimensions. This is shown in full for 1- and 2-dimensional strain below:


One - dimensional Lagrangian strain. The length is the only strain component, and thus L is measured along the only coordinate axis, thus L = x



Strain in two dimensions. Above are the two normal strains along the x and y axes, where each strain component can be seen as Lagrangian strain along one main axis.  Below are the two shear strain components, movement of the borders relative to each other. Here there are two strain components, characterised by the tangent to the shear angle alpha.

Thus in two dimensions, the strain tensor has four components, two normal:


And two shear components:

From the figure, it is also evident that:

The whole strain tensor can be written as a matrix:



In three dimensions there are three normal and six shear strain components:


as shown below:

Strain in three dimensions. Only the three strain components along the x axis  (one normal, two shear) are shown, but the y and z strains will be exactly the same and can be imagined by rotating the x images. The three strain components shown are:


Coordinate systems:

Strain in three dimensions has so far been discussed in terms of a traditional Cartesian xyz coordinate system, for general purposes. But three dimensions can be defined by any coordinate system that defines a point in space by three coordinates. A coordinate system can be defined in relation to the ultrasound imaging plane: Axial (depth - i.e. along the ultrasound beams), lateral (In-plane angle or distance - i.e. across the beam) and elevation (out of plane distance or angle).


Ultrasound coordinates: The coordinate system gives the position relative to the ultrasound imaging plane. The three directions axial, lateral and elevation can be defined as  orthogonal coordinates (green), but one that is relative to the position of the imaging plane instead of fixed in space. However, both lateral and elevation can be defined by the angle instead, i.e. a coordinate system of one distance and two angles.

In relation to the left ventricle, a local coordinate system is most used; longitudinal, transmural and circumferential. At any given position, those coordinates are orthogonal, but the direction of the circumferential and transmural axis changes as one moves around the ventricle (so does longitudinal if the direction is defined along the wall instead of along the main axis of the ventricle):


(Figure courtesy of A Heimdal)


Incompressibility

That an object is incompressible means that the volume (not mass!) remains constant during deformation as shown below:

Incompressibility. The figure shows simultaneous longitudinal shortening or compression, as well as thickening or expansion in the two transverse directions. If the cylinder is incompressible, the sum of the longitudinal and the two transverse strains will be zero.  The cube increases its length along the x axis (positive strain), while the x and y lengths decreases (negative strain), in such a way as to keep the volume constant:



The volume of the cube is:
V = x × y × z
After deformation, the volume is:

As the volume remains constant during deformation, this means that:

hence:

and:

and thus:



As evident from this formula, there can be both shear strains and normal strains in an object, i.e. 3D deformation, as long as they even out.


In the heart, the incompressibility of the myocardium means that the volume of the myocardium remains constant. As the outer contour remains constant during systole as discussed in long axis function, this means that as the ventricle shortens, the wall has to thicken inwards.
As the outer contour remains constant during systole as discussed in long axis function, this means that as the ventricle shortens, the wall has to thicken inwards. When the left ventricle shortens in systole, the total volume is reduced by the volume of the cylinder  shown in grey.
 The left ventricle seen as a half ellipsoid with apex wall thickness half the wall thickness at the base. For the left ventricle, the short axis cross section may be assumed circular, and the short axis radii thus equal (r = d/2), while the long axis radius is equal to the length of the half ellipsoid, l. Wall thickness at the base; w is the difference between outer and inner radius (or two wall thicknesses is the difference between outer and inner diameter): do = di + 2W (or inner diameter di = do - 2W). 
Wall thickness decreases toward the apex, for simplified calculation (to avoid recurion) it is assumed to be zero. li = lo = l

The volume of a general ellipsoid is 4/3 * pi * abc, where a, b and c are the three radii of the ellipsoid.  In this case with a circular cross section; a = b = r = d/2 while c is the long radius and c = l. Thus the volume of a half ellipsoid in this case is 2/3* pi * (d/2)2 * l.
Total or outer volume Vt = 2/3* pi * (do/2)2 * l and cavity (inner) volume Vc = 2/3* pi * (di/2)2 * l. Wall volume is Vw = Vt - Vi
 
Midwall diameter is dm = di + W and midwall circumference is cm  = dm * pi.


Longitudinal strain;  l, thus systolic length; ls = ld + ld * l



 Diastole
 Systole
Total volume Vtd = 2/3* pi * (do/2)2 * ld Vts = 2/3* pi * (do/2)2 * (ld + ld *l)   With an invariant outer diameter, this would equal Vtd + Vtd * l = Vtd (ld+  l)
Cavity volume diastole
 Vcd = 2/3* pi * (did/2)2 * ld  
Wall volume
 Vw = Vwd = Vtd - Vcd Vws = Vwd = Vw due to incompressibility
Cavity volume systole
Vcs = Vts -Vw
Stroke volume
SV = Vcd - Vcs 
EF
EF = SV / Vcd  With an invariant outer diameter, this would equal (-Vtd * l) / (Vtd - Vw)
Wall thickness
 Wd = 1/2 (do - did)
 Ws = Wd + Wd * WT
Fractional shortening
FS = (did - dis) / dis

Thus it is evident that cavity volume is a function of LV length, outer diameter and wall thickness.
Under the incompressibility condition and the assumption of an invariate outer diameter, longitudinal shortening (strain) is the sole determinant of wall thickening, (transmural strain) and hence, midwall circumferential strain, stroke volume and EF, as well as fractional shortening. Wall thickening is given by wall shortening in order to keep wall volume constant, but does not enter into the equations for stroke volume.

However, as the circular fibres contract, the outer diameter may change. This will contribute to both wall thickening as well as well as inner diameter reduction. In this case also to stroke volume and EF. 


In  this model, a left ventricle with diastolic length 9.5 cm, outer diameter of 7 cm and wall thickness of 10 mm (reasonable findings in accordance with the findings in the HUNT study (153)), a 22% longitudinal strain will result in:


Invariant outer diameter
 5% systolic reduction in outer diameter
Stroke volume
 44 ml
 72 ml
EF
 33%
 58%
Wall thickening
 27%
 49%
Midwall shortening
 4%
 14%
Fractional (inner diameter) shortening
 9%
 27%

So it seems that the eggshell concept is slightly inaccurate, and a ca 5% reduction in outer diameter, close to that found in clinical studies (158), also results in more normal functional measures.





Strain in three dimensions versus EF








Angle deviation of velocity, motion and strain rate

(Se also problems and pitfalls, insonation angle deviation.) A more advanced treatment is given here.

It is well known that velocity measurement is dependent on the angle between the ultrasound beam and the velocity direction (vector) –insonation angle, and that the measured velocity vm is reduced, compared to the true velocity v by the cosine of the insonation angle.



Left: distance along the axis x, imagined along the ultrasound beam y. The angle effect is:


This means that motion along an M-mode line, actually will increase by the cosine of the angle deviation, as opposed to the velocity. This is explained to the right. The angle deviation gives an increase in wavelength, similar to the increase in distance:

Frequency of ultrasound is the inverse of wavelength, so:

And thus:


This results in a 6% underestimation at an angle of 20º, 13% at 30º, 29% at 45º and zero velocity at 90º.

Basically, this is an oversimplification. Points on the outer contour of the myocardium will move longitudinally, while endocardial points will move inwards as well. Thus, the endocardial points has a transverse velocity component as well. This means that even though each point has only one velocity vector, the angle between the vector and the longitudinal direction of the wall (and hence the insonation angle as well), will vary across the wall, not only with alignment of the ultrasound beam with the wall, but with position across the wall as well:



 For strain rate and strain measurement, the angle dependency is somewhat more complex. In the longitudinal direction, the longitudinal velocity gradient is reduced by the cosine of the insonation angle, as for velocities. But in an incompressible object, there is simultaneous strain in the transverse direction, in order to keep the volume constant, and the two strain tensors are opposite and will detract from each other (1, 2, 7). This results in further reduction in the measured strain and strain rate (7). The cosine formula above, may explain some of the angle relations in a simplified way:



Longitudinal velocity and motion are vl and xl, respectively. Transverse velocity and motion are vt and xt.  Along the ultrasound beam, the longitudinal and transverse velocity components are:

 and

the total velocity is the sum of the two components:

The total distance is the diagonal in the rectangle:

The formula is valid for velocity differences, as can be shown by setting v2 = 0. Then:






It must be emphasised that this argument is simplified, to illustrate the principles, and not an algorithm for actual angle correction. Firstly, in the heart there are three dimensions, not two: longitudinal, transmural and circumferential. In addition, the interrelation between all three components is of importance. For a more complete mathematical analysis, see (1). The angle problem is also specially analysed in: Andreas Heimdal. Angle dependency of systolic strain measurements using strain rate imaging.


In addition there are shear strain tensors and finally there are angle related artefacts that are not subject to mathematical analysis, which may be even more important as discussed here.


Back to website index
Back to section index

References:


Editor: Depertment head, Contact address: isb-post@medisin.ntnu.no, Updated: XXX