Japanese Journal of Applied
Physics

You may also like

REGULAR PAPER

High-speed measurement of two-dimensional
displacement of myocardium using element RF
data of ultrasonic probe
To cite this article: Kaisei Hara et al 2023 Jpn. J. Appl. Phys. 62 SJ1040

View the article online for updates and enhancements.

- Power cavitation-guided blood-brain
barrier opening with focused ultrasoundand microbubbles
M T Burgess, I Apostolakis and E E
Konofagou- Locally optimized correlation-guided
Bayesian adaptive regularization for
ultrasound strain imaging
Rashid Al Mukaddim, Nirvedh H Meshram
and Tomy Varghese- Bone surface enhancement in ultrasoundimages using a new Doppler-based
acquisition/processing method
Xu Yang, Songyuan Tang, Ennio Tasciotti
et al.

This content was downloaded from IP address 130.34.32.223 on 12/04/2023 at 04:44

REGULAR PAPER

Japanese Journal of Applied Physics 62, SJ1040 (2023)
https://doi.org/10.35848/1347-4065/acc112

High-speed measurement of two-dimensional displacement of myocardium using
element RF data of ultrasonic probe
Kaisei Hara1, Shohei Mori2*, Mototaka Arakawa1,2, and Hiroshi Kanai1,2
1Graduate School of Biomedical Engineering, Tohoku University, Sendai, Miyagi 980-8579, Japan
Graduate School of Engineering, Tohoku University, Sendai, Miyagi 980-8579, Japan2*

E-mail: mori@tohoku.ac.jp

Received November 15, 2022; revised February 16, 2023; accepted March 2, 2023; published online March 27, 2023
Improving the accuracy of heart wall motion measurement is essential to realise better cardiac function evaluation. This paper proposed a twodimensional (2D) displacement estimation method with a high temporal resolution using the 2D complex cross-correlation of element RF signals of
an ultrasonic probe between frames returned from the target scatterers. The application of the proposed method to the phantom displacement
confirmed its principle. The estimated 2D displacement of the phantom was consistent with the set displacement. Subsequently, the method wasapplied to two healthy subjects to measure the 2D displacement of the interventricular septum during one cardiac cycle. Consequently, duringsystole and diastole, the movement of the myocardium was measured, and the results were validated.
© 2023 The Japan Society of Applied Physics

1. Introduction

Echocardiography is useful for evaluating cardiac function
owing to its repeatability, non-invasive diagnostic capability,and high temporal resolution. Consequently, various medical
applications have been studied to evaluate cardiac function,
including measuring heart wall strain rates1–7) and propagation of response of the myocardium to contraction.8–10)
However, further enhancement of the accuracy of heart
wall motion measurements is essential to realise better
cardiac function evaluations.
The speckle tracking method was applied to measure thetwo-dimensional (2D) displacement between frames of theheart wall using the correlation window on the acquired Bmode image to track changes in the speckle position within
the correlation window between consecutive frames.11–17) By
increasing the number of transmissions of the beam direction
per frame, smaller sampling intervals in the lateral direction
were achieved; however, the temporal resolution was reduced. In clinical practice, heart wall displacement has beenmeasured with a temporal resolution in the slow range of 40–
80 Hz for normal heart rates.18–20)
However, in vivo measurements of the instantaneous displacement (velocity) along the ultrasonic beam of the heart wall
revealed the existence of higher frequency components.
Moreover, for these measurements, a higher frame rate(>400 Hz) is necessary to improve the accuracy of the velocity
measurement.10) Consequently, 2D displacement measurements
require a similar high frame rate to maintain the accuracy.
To measure the 2D displacement of the heart wall with ahigh temporal resolution, a method for forming an arbitrary
number of received beams per transmission of unfocused
waves was proposed.21–25) As unfocused wave transmission
realises a spatially wider range of measurements, thistechnique signiﬁcantly improves the temporal resolution
and achieves a high density of scan lines in the lateral
direction.26,27) However, signals from other scatterers thatmove independently are also added and spatially averaged,
which contaminates the amplitude of the receiving beam andaffects the measurement of the 2D dimensional displacement.
The received time distribution of the RF waves of theprobe elements scattered from the target changes according to

the 2D displacement of the target. A method for measuring
the 2D displacement without beam scanning was proposed bymeasuring the change in the received time distribution of thescattered waves.28,29) However, these previous studies did notconsider the actual distribution of the received time on probeelements. The received time distribution was approximatedby using a linear function determined by the element
position. These studies aimed to measure the 2D velocity
vector of blood ﬂow; however, the measurement of the 2D
displacement of the heart wall was not fully considered.
The difference in the received time distribution is small
because the 2D displacement of the heart wall between
frames is smaller than that for blood ﬂow. Thus, considering
the actual distribution of the received time on probe elements
is an effective strategy. Based on this principle, a method toestimate the 2D displacement by evaluating the root meansquared difference of the received time distribution from thegeometrically calculated received time distribution was
proposed.30) However, the condition with multiple scatterers
was not considered.
Thus, this study proposes a novel 2D displacement
estimation method with a high temporal resolution. The 2D
displacement was estimated using the 2D complex crosscorrelation of RF signals of probe elements between frames,
wherein RF signals were arranged by the received time
distribution of scattered waves on the probe. The applicationof the proposed method to the phantom displacement
conﬁrmed its principle from the estimate of the 2D displacement of the phantom. Thereafter, the proposed method wasapplied to two healthy subjects and the 2D displacement of
the interventricular septum (IVS) during one cardiac cycle
was successfully measured.
2. Principles
2.1. 2D displacement estimation method using
received time distribution of scattered waves on
elements of ultrasonic probe

As illustrated in Fig. 1, consider the case wherein the target
scatterer was set at (xs, zs) at the nth frame and was displacedby (Δx, Δz) at the (n + 1)th frame, where Δx and Δz denotethe instantaneous displacements in the lateral and axial
directions, respectively, between two consecutive frames.

SJ1040-1© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

Fig. 3. (Color online) Schematic of the water tank experimental systemwith the agar phantom.

zs2 + (xs - xk )2
zs
+
,
c0
c0
K
K
⎛for k = - , ¼, 0, ¼, ⎞2
2⎠⎝

T (x k ; x s , z s ) =Fig. 1. (Color online) Time distribution {T (xk ; xs , zs )} of scattered waves
received by the probe elements after the plane wave is transmitted. The blue
arrows show the propagation path of the plane wave from the probe to thetarget scatterer, and the red arrows show the propagation path from the target
scatterer at the position (xs , zs ) to each element at the n th frame. The target
scatterer is displaced by (Dx, Dz ) at the (n + 1)th frame.

(a)

(b)

(1 )

where (K + 1) is the number of probe elements.
Figure 2(a) illustrates the time distribution {T (xk ; xs, zs)}of the RF signals on the probe elements. The region withlength ΔT between the two red lines shows the duration timeof the scattered waves from the target, where the region starts
at the ﬁrst arrival time T (xk ; xs, zs). Figure 2(b) shows the{T (xk ; xs + Dx, zs + Dz )} scattered from the target at(xs + Dx, zs + Dz ) at the (n + 1)th frame. Owing to theinstantaneous 2D displacements (Dx, Dz ) of the measurement target, the time distribution of the scattered waves
changed from red to black lines.
For a sufﬁciently high frame rate, the change in the
waveforms of the waves scattered from the target between
the two frames can be assumed to be small. Thus, theassumed instantaneous 2D displacements (Dx, Dz ) of thetarget between the frames can be evaluated by the following
2D complex cross-correlation Cn (Dx, Dz; xs, zs) using thecorrelation window with length DT:
Cn (Dx , Dz ; xs , zs)
1Fig. 2. (Color online) Principle of the proposed method to estimate 2D
displacement (Dx, Dz ) of the target from n th frame to the (n + 1)th frame.
(a) Time distribution {T (xk ; xs , zs )} of RF signals of probe elements and theregion with the length DT between the two red lines, where DT corresponds
to the duration time of the scattered waves. (b) Time distribution{T (xk ; xs + Dx, zs + Dz )} scattered from the target at (xs + Dx, zs + Dz )at the (n + 1)th frame (black line).

Let the centre element position at x0 = 0 among (K + 1)
elements be the origin. After the plane wave is transmitted
from all elements of the probe at time 0, it is scattered by thetarget at positions (xs, zs), and then the RF signal is received
by each element of the probe. The ﬁrst arrival timeT (xk ; xs, zs) of the scattered wave received by the kthelement at xk is expressed as the round-trip propagation
distance divided by the speed c0 of sound in the medium asfollows:

=A

K /2

DT

å å yn*,k (T (xk ; xs, zs) + t )k =-K / 2 t = 0(2 )·yn + 1, k (T (xk ; xs + Dx , zs + Dz) + t ) ,where yn, k (t ) is the complex waveform obtained via theapplication of quadrature demodulation to the waveform
acquired by the k th element in the n th frame, * denotes thecomplex conjugate, and A is the multiple of the root meansquares of RF signals in the correlation window at the n th and(n + 1)th frames expressed as
A=

K /2

DT

å å ∣ yn,k (T (xk ; xs, zs) + t )∣2

k =-K / 2 t = 0
K /2·

SJ1040-2DTå å ∣yn + 1,k (T (xk ; xs + Dx , zs

(3 )+ Dz) + t )∣2 .

k =-K / 2 t = 0© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

Fig. 4. (Color online) B-mode image of a parasternal left ventricular longitudinal view obtained by the ultrasound diagnosis apparatus for subject A. RV,
IVS, LV, ECG, and PCG are the abbreviations of the right ventricle, the interventricular septum, the left ventricle, electrocardiogram, and phonocardiogram(heart sound), respectively.

water
agar
phantom

Fig. 5. B-mode image of the phantom in Fig. 3 with a width of the lateralview of 50° and the positions of ﬁve targets s1, s2, s3, s4, and s5 set in the agar
phantom.

The instantaneous 2D displacement of the target is deterx , Dz ) thatmined by ﬁnding the assumed displacements (Dmaximise the absolute value of Cn (Dx, Dz; xs, zs) as
x , Dz) = arg
(Dmax

(Dx , Dz ) Î R∣ Cn (Dx , Dz ; xs , zs)∣ ,(4 )where R denotes the range for search in axial and lateral
directions.
When the non-target scatterer is set at (xs0 , zs0 ) at the nthframe, which is then displaced by (Dxs0 , Dzs0 ) at the (n + 1)
th frame, the scattered waves from the target scatterer at
positions (xs, zs) and those from the non-target scatterer at
positions (xs0 , zs0 ) are overlaid. However, the time difference
distribution of the scattered waves from the target scatterer,{T (xk ; xs + Dx, zs + Dz )} - {T (xk ; xs, zs)}, differs fromthat of the scattered waves from the non-target scatterer,{T (xk ; xs0 + Dxs0 , zs0 + Dzs0 )} - {T (xk ; xs0 , zs0 )}. Thus,
it is expected that the absolute values ofCn (Dx, Dz; xs, zs) of Eq. (2) exhibits the highest value atthe position corresponding to the 2D displacement (Dx, Dz )for the target.
2.2. Phantom experiment
To verify the accuracy of the proposed method, it was used tomeasure 2D displacement of a phantom in a water tankexperiment, where the phantom was created by mixing agar
(concentration of approximately 10%) with water. Although
the agar phantom did not include a scatterer, we conﬁrmedthat scattered waves from inside the agar phantom werereceived based on its inhomogeneity. The phantom with 26
and 35 mm in the axial and lateral directions, respectively as
shown in Fig. 3, was measured. The phantom was set in awater tank, and the probe was moved in the x–z stage. The
displacement of the probe was set as Δz = −0.2 mm in theaxial direction and Δx = 1 mm in the lateral direction, whichwere determined based on the maximum instantaneous 2D
displacement of the human myocardium at the IVS between
the consecutive frames during the rapid ﬁlling phase at aframe rate of 200 Hz.
The element RF signals were acquired after the transmission of the single-plane wave using an ultrasound diagnostic
apparatus (Prosound SSD-α10, Hitachi Aloka, Japan). The
number of elements in the sector probe was K = 95, with anelement spacing of 0.2 mm (UST-52101N, Aloka). The
centre and sampling frequencies of the transmitted wave
were 3.75 and 30 MHz, respectively. Since the watertemperature was 26.7 °C, the speed of sound was calculated
to be 1501 m s−1.31)
To calculating the 2D complex cross-correlation
Cn (Dx, Dz; xs, zs), the length DT of the correlation window
was set to 0.8 μs (1.2 mm). To determine (Dx, Dz ) in Eq. (4),the range R of ±1.2 mm was searched in the axial and lateral
directions.

SJ1040-3© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)(a)

K. Hara et al.

(b)

Fig. 6. (Color online) (a) Amplitudes of the element signals acquired before displacing the probe in Fig. 3. The upper red line shows the received time
distribution of scattered waves, {T (xk ; xs1, zs1)}, for the target s1 at (xs1, zs1) = (-5 mm, 56 mm). The lower red line shows the end time of the correlation
window {T (xk ; xs1, zs1) + DT}. (b) Amplitudes of the element signals acquired after displacing the probe’s position by (Dx, Dz ) = (1 mm, -0.2 mm). The
x, zs + Dz )}red dashed lines: the distributions {T (xk ; xs1, zs1)} and {T (xk ; xs1, zs1) + DT} in Fig. 6(a). The black lines show the distributions {T (xk ; xs1 + D
1
x, zs + Dz ) + DT}.
and {T (xk ; xs1 + D
1

RV
IVS

LV
Fig. 7. (Color online) Magnitude distribution of the 2D complex crosscorrelation, ∣ Cn (Dx, Dz; xs1, zs1)∣, obtained from element signals in
Fig. 6(a) to those in Fig. 6(b).

Fig. 9. For subject A, the B-mode image of a parasternal left ventricularlongitudinal view with a width of the lateral view of 50° at the initial frame of
the analysis (−100 ms from TR for the ﬁrst heartbeat). The position(xs , zs ) = (2.07 mm, 39.57 mm), where the 2D displacement of the IVS
was measured, is shown by a white cross mark.

2.3. In vivo experimentFig. 8. (Color online) For each of the ﬁve targets s1, s2, s3, s4, and s5, theposition was set at the white cross mark before the displacement of the probe,
as shown in Fig. 5. After the displacement of the probe, the true position of
the 2D displacement (white circle) and the estimate (red cross mark) are
shown for each target.

The proposed method was applied to the measurement of the
2D displacement of the IVS during one cardiac cycle and theresults for each cardiac phase were discussed. The subjects
were two healthy males in their 20s. Figure 4 shows a Bmode image of the parasternal left ventricular longitudinalview of subject A. The element RF signals from the IVS wereacquired each time a plane wave was transmitted immediately below the probe. The speed of sound was set to1540 m s−1, and the frame rate was set to 1205 Hz, that is,
the frame interval DTs ¢ = 830 ms. Although the theoreticallyachievable frame interval was DTs = 166 ms for the employed ultrasound diagnosis equipment, DTs ¢ = 830 ms was
set in the experiment owing to the transmission conditions forthe equipment.
For calculating Cn (Dx, Dz; xs, zs), the length DT of thecorrelation window was set to 0.78 ms (1.2 mm). Further, to

SJ1040-4© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

Fig. 10. (Color online) (I) Maximum absolute value of the 2D complex
cross-correlation ∣ Cn (Dx, Dz; xs , zs )∣ for each frame n, (II) estimated xˆ (n )of the lateral displacement, (III) estimated zˆ (n ) of the axial displacement,(IV) ECG, and (V) PCG, during one cardiac cycle for subject A. For theseven typical times (a)–(g), element RF signals and the distribution of the 2D
complex cross-correlation are shown in Fig. 11.

determine (Dx, Dz ) in Eq. (4), a range R of ±0.4 mm wasset in the axial and lateral directions. These measurements
were approved by the Ethics Committee of the Graduate
School of Engineering, Tohoku University, and all theparticipants provided consent before participating in thepresent study.
3. Results and discussion
3.1. Phantom experiment

Figure 5 shows the B-mode image of the phantom in Fig. 3
with a lateral view width of 50°, and a solid closed line
surrounding the phantom. Five measurement positions, s1, s2,
s3, s4, and s5 were set in the phantom, as shown in Fig. 5.
Figure 6(a) shows the amplitudes of the element signals
acquired before displacing the probe (Fig. 3). For target s1 at(xs1, zs1) = (-5 mm, 56 mm), the received time distributionof scattered waves, {T (xk ; xs1, zs1)}, was determined using
Eq. (1), as indicated by the red upper line. The end time
distribution of the correlation window, {T (xk ; xs1, zs1) + DT },is also indicated by the lower red line.
Because the phantom surface was not ﬂat, as shown in
Fig. 3, and the phantom was closest to the probe at

xk = −5 mm, the lateral measurement position x si was set
to −5 mm for all targets (i = 1, …, 5). For targets s2, s3, s4,
and s5, the axial positions {z si} were set apart by 5 mm in therange of 61– 76 mm in depth (Fig. 5).
Figure 6(b) shows the amplitudes of the element signals
acquired after displacing the probe position by (Dx, Dz ) =(1 mm, -0.2 mm). The distributions {T (xk ; xs1, zs1)} and {T (xk ; xs1, zs1) + DT } in Fig. 6(a) are indicated by red dashedx , Dz ), the distributionslines. For the assumed positions (Dz ) + DT} are{T (xk ; xs1 + Dx, zs1 + Dz )} and {T (xk ; xs1 + Dx, zs1 + D
represented by black lines, where (Dx , Dz ) are determinedby the position with the maximum correlation value (Fig. 7).
A comparison of Figs. 6(a) and 6(b) revealed that thereceived time distribution of the scattered waves shifted
almost parallel from the positions of the red lines to those of
the black lines according to the 2D displacement of(Dx, Dz ) = (1 mm, -0.2 mm).
Figure 7 shows the magnitude distribution of the 2D
complex cross-correlation, ∣ Cn (Dx, Dz; xs, zs)∣, of Eq. (2),from the element RF signals in Fig. 6(a) to those in Fig. 6(b).
From the position with the maximum correlation value of∣ Cn (Dx, Dz; xs, zs)∣ in Fig. 7, the instantaneous 2D displax , Dz ) =
cement of the target s1 was estimated to be (D(-0.98 mm, 0.23 mm), which was consistent with the set
displacement (-Dx, -Dz ) =(-1 mm, +0.2 mm) for thephantom as described above.
Figure 8 shows the 2D displacement estimate for each of
the targets s1, s2, s3, s4, and s5 in Fig. 5, which is indicated by
the red cross mark. The white circle indicates the true
position after the 2D displacement of each target. The
mean and standard deviation of the 2D displacement of theseﬁve targets in Fig. 5 for the lateral and axial directions were
x , Dz ) = (-1.004  0.023 mm, 0.234  0.008 mm).
(D
As shown in Fig. 8, the 2D displacement estimated for eachtarget matched well with the set displacement(-Dx, -Dz ) = (-1 mm, +0.2 mm).
3.2. In vivo experiment
Figure 9 shows the B-mode image of a parasternal left
ventricular longitudinal view for subject A with a lateral view
width of 50° at the initial frame of the analysis (−100 ms
from the time TR of the R-wave for the ﬁrst heartbeat). The
target was set in the IVS at position (xs, zs) =(2.07 mm, 39.57 mm) , where its 2D displacement wasmeasured below. The target was set by selecting the brightest
position in the B-mode image.
Figure 10 shows the (I) maximum value of the∣ Cn (Dx, Dz; xs, zs)∣ for each frame n, (II) estimated xˆ (n )of the lateral displacement, (III) estimated zˆ (n ) of the axial
displacement, (IV) electrocardiogram (ECG), and (V) phonocardiogram (PCG) during one cardiac cycle for subject A.
As shown in Figs. 10(II) and 10(III), each displacement
waveform was reset to 0 at TR for the ﬁrst heartbeat (time 0).
Seven typical times (a)–(g) were chosen (Fig. 10):(a) initial frame, (b) beginning of the isovolumetriccontraction phase (TR for the ﬁrst heartbeat), (c) beginningof the ejection phase (the time of the ﬁrst heart sound),(d) beginning of the isovolumetric relaxation phase (thetime of the second heart sound), (e) the rapid ﬁlling phase,(f) beginning of the atrial systole (the time of P-wave
for the second heartbeat), and (g) TR for the second
heartbeat.

SJ1040-5© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

element position

element position

Fig. 11. (Color online) (left) Amplitudes of the element signals with {T (xk ; xs , zs )} and {T (xk ; xs , zs ) + DT} for two consecutive frames and (right)∣ Cn (Dx, Dz; xs , zs )∣, at the seven typical times (a) to (g) for subject A. (a) Initial frame, (b) beginning of the isovolumetric contraction phase (TR for the ﬁrst
heartbeat), (c) beginning of the ejection phase, (d) beginning of the isovolumetric relaxation phase, (e) the rapid ﬁlling phase, (f) beginning of the atrial systole,and (g) TR for the second heartbeat.

As shown in Figs. 10(II) and 10(III), the heart wall
rapidly moved toward the apical side and the left ventricle
from (b) to (c), which corresponds to the movement of the
myocardium during systole. Following the movement of

the heart wall toward the basal side and the right ventricle
from (d) to (f), it moved toward the left ventricle, whichcorresponded to the movement of the myocardium during
diastole.

SJ1040-6© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.
(b)

(a)

Fig. 12. (top) For subject A, the B-mode image of a parasternal left ventricular longitudinal view and (bottom) B-mode image magniﬁed near themeasurement position (xs , zs ). (a) One frame before when max (Dx, Dz) Î R ∣ Cn (Dx, Dz; xs , zs )∣ was lowest (381.8 ms from TR for the ﬁrst heartbeat),(xs , zs ) = (-4.66 mm, 36.87 mm). (b) Time when max (Dx, Dz) Î R ∣ Cn (Dx, Dz; xs , zs )∣ was lowest (382.63 ms from TR for the ﬁrst heartbeat),(xs , zs ) = (-4.61 mm, 36.87 mm).

Figure 11 shows the amplitudes of the element signals with{T (xk ; xs, zs)} and {T (xk ; xs, zs) + DT} for two consecutive
frames and ∣ Cn (Dx, Dz; xs, zs)∣ at seven times (a)–(g).
From the initial frame to the beginning of the ejection phase,as shown in Figs. 11(a)–11(c), the correlation distribution
was unimodal. However, from the isovolumetric relaxation
phase to the rapid ﬁlling phase, as shown in Figs. 11(d) and11(e), the amplitude of the element RF signals decreased, andthe correlation distribution then became lower and obscured.
For atrial systole, as shown in Figs. 11(f) and 11(g), both theamplitude of the element RF signal and the correlation
increased again.
Figure 12(a) shows the B-mode image of a parasternal left
ventricular longitudinal view for subject A at one frame
before when max(Dx, Dz) Î R ∣ Cn (Dx, Dz; xs, zs)∣ was lowest(381.8 ms from TR for the ﬁrst heartbeat) during the rapid
ﬁlling period between (d) and (e) in Fig. 10 and its magniﬁed
image around the measurement position (xs, zs) =(-4.66 mm, 36.87 mm). Figure 12(b) shows the B-mode
image of the same view when max(Dx, Dz) Î R ∣Cn (Dx, Dz; xs, zs)∣
was lowest (382.63 ms from TR for the ﬁrst heartbeat) and its

magniﬁed image around the position (xs, zs) =(-4.61 mm, 36.87 mm). These B-mode images were obtained by normalising the envelope amplitude of the received
beams to the same value as in Fig. 9. Comparisons of themagniﬁed images in Figs. 12(a) and 12(b) showed that the
speckle pattern near the position (xs, zs) in the IVS changed
during DTs ¢ = 830 ms between two consecutive frames,
which would be caused by the rapid thinning of the heartwall. Thus, the decrease in the correlation distribution shownin Figs. 11(d) and 11(e) was caused by a change in waveforms owing to the change in thickness of the heart wall.
Figure 13 shows the maximum values of∣ Cn (Dx, Dz; xs, zs)∣ for the same data set at the initialframe of Fig. 10(a) for various frame rates (the time intervalfrom the initial frame). The maximum values signiﬁcantly
decreased for a frame rate of less than 100 Hz owing to thechange in waveforms between frames with lower time
resolution. However, for a high temporal resolution(=1205 Hz) as achieved in the proposed method, thecorrelation increased to estimate the instantaneous 2D
displacement.

SJ1040-7© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

Fig. 13. (Color online) Maximum values of ∣ Cn (Dx, Dz; xs , zs )∣ for thesame data set at the initial frame of Fig. 10(a) for various frame rates (thetime interval from the initial frame).

RV

IVS
LV

Fig. 14. For subject B, the B-mode image of a parasternal left ventricularlongitudinal view at the initial frame of the analysis (−100 ms from TR forthe ﬁrst heartbeat). The position (xs , zs ) = (-6.02 mm, 43.2 mm), where
the 2D displacement of the IVS was measured, is shown by a white crossmark.

For subject B, Fig. 14 shows the B-mode image of aparasternal left ventricular longitudinal view at the initialframe of the analysis (−100 ms from TR for the ﬁrst heartbeat). The target was set in the IVS at position(xs, zs) = (-6.02 mm, 43.2 mm), where the target was set
by selecting the brightest in the B-mode image.
Figure 15 shows the (I) maximum value of the∣ Cn (Dx, Dz; xs, zs)∣ for each frame n, (II) estimates xˆ (n )of the lateral displacement, (III) estimates zˆ (n ) of the axial
displacement, (IV) ECG, and (V) PCG, during one cardiac
cycle for subject B. As shown in Figs. 15(II) and 15(III), thetarget moved toward the apical side and the left ventricle
from (c) to (d) during the systole, and it moved toward thebasal side and the right ventricle from (d) to (g) during thediastole, similar to subject A as shown in Figs. 10(II) and10(III). Further, as shown in Figs. 10(II), 10(III), 15(II), and15(III), the maximum displacements in the lateral direction
and the axial direction for subjects A and B were(
x max , zmax ) = ( - 9.64 mm, - 3.2 mm) and

Fig. 15. (Color online) (I) Maximum absolute value of the 2D complex
cross-correlation ∣ Cn (Dx, Dz; xs , zs )∣, (II) estimates xˆ (n ) of the lateraldisplacement, (III) estimates zˆ (n ) of the axial displacement, (IV) ECG, and(V)PCG, during one cardiac cycle for subject B. For the seven typical times(a)–(g), element RF signals and the distribution of the 2D complex crosscorrelation are shown in Fig. 16.

(-10.4 mm, -4.76 mm), respectively, which were in thesame order as the previous study (-8 mm, -3 mm).32,33)
Figure 16 shows the amplitudes of the element signals with{T (xk ; xs, zs)} and {T (xk ; xs, zs) + DT} for two consecutive
frames and ∣ Cn (Dx, Dz; xs, zs)∣ at seven typical times (a) to(g), as shown in Fig. 15 for subject B. From the isovolumetricrelaxation phase to the rapid ﬁlling phase, as shown in
Figs. 16(d) and 16(e), the amplitude of the element RF
signals decreased. Subsequently, the correlation distribution
reduced and obscured, similar to that observed in subject A,
as shown in Figs. 11(d) and 11(e). For the atrial systole, as
shown in Figs. 16(f) and 16(g), both the amplitude of theelement RF signal and the correlation increased again, similar
to that observed in subject A, as shown in Figs. 11(f) and11(g).
3.3. Comparison between the proposed and speckle
tracking methods
To demonstrate the superiority of our method, we compared
the proposed method with previous studies using the speckle
tracking method and received time distribution. First, we
compared the proposed method with a conventional speckle
tracking method.

SJ1040-8© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Fig. 16. (Color online) (left) Amplitudes of the element signals with {T (xk ; xs , zs )} and {T (xk ; xs , zs ) + DT} for two consecutive frames and (right)∣ Cn (Dx, Dz; xs , zs )∣, at the seven typical times (a) to (g) for subject B. (a) Initial frame, (b) beginning of the isovolumetric contraction phase (the time TR of
the R-wave for the ﬁrst heartbeat), (c) beginning of the ejection phase, (d) the beginning of the isovolumetric relaxation phase, (e) the rapid ﬁlling phase, (f)beginning of the atrial systole, and (g) TR for the second heartbeat.

x (n ) of the instantaFigure 17 shows the (I) estimated D
z (n ) of theneous lateral displacement, (II) estimated D
instantaneous axial displacement, (III) ECG, and (IV) PCG,
during one cardiac cycle for subject A.

For the conventional speckle tracking method, uponassuming 112 beams being employed to obtain the 2D
envelope signals with a width of the lateral view of 42°, aperiod of 112 ´ DTs = 18.59 ms was necessary to acquire

SJ1040-9© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.
(a)

(5)
(4)

estimated position

(3)
(1)

(2)

estimated position(b)

instantaneous axial displacement(3)x (n ) of the instantaneous lateralFig. 17. (Color online) (I) Estimated D
z (n ) of the instantaneous axial displacement,
displacement, (II) estimated D(III) ECG, and (IV) PCG, during one cardiac cycle for subject A.

the 2D signal when the transmitting interval DTs = 166 ms.
For two consecutive frames, twice of 18.59 ms, that is,37.18 ms was necessary, which provides the time resolutionof the conventional speckle tracking method and it wasapproximately 44 times longer than the frame interval
employed in the present experiment, DTs¢ = 830 ms.
Figure 18(a) shows the transition in the estimated positions{(xs (n ), zs (n))} of the target for 45 frames fromTR - 18.26 ms to TR + 18.26 ms in Fig. 17. The ﬁve typicaltimes (1)–(5) are shown by the red circles, that is, (1)TR - 18.26 ms, (2) TR - 9.13 ms, (3) TR, (4) TR + 9.13 ms,
and (5) TR + 18.26 ms.
For this period, the transition in the estimated 2Dx (n ), Dz (n ))} in Fig. 17
instantaneous displacements {(Dis shown in Fig. 18(b), and the estimates at the ﬁve typicaltimes are shown by red circles. As shown in Fig. 18(b), the2D instantaneous displacement rapidly changed in short

(4)
(5)
(2)

(1)
instantaneous lateral displacement

Fig. 18. (Color online) (a) Transition in the estimated positions{(xs (n ), zs (n ))} of the target for 45 frames from TR - 18.26 ms to
TR + 18.26 ms. The ﬁve typical times (1) to (5) were shown by the redcircles, that is, (1) TR - 18.26 ms, (2) TR - 9.13 ms, (3) TR, (4)TR + 9.13 ms, and (5) TR + 18.26 ms. (b) Transition in the estimated 2Dx (n ), Dz (n ))} in Fig. 17. The estimates atinstantaneous displacements {(Dthe ﬁve typical times were shown by the red circles.

3.4. Comparison of the proposed method with thoseof previous studies using received time distribution

Next, we compared our method with previous studies28,29)
using the received time distribution of scattered waves on theelements in the ultrasonic probe. The time difference
distribution DT (xk ; xs, zs, Dx, Dz ) between T (xk ; xs, zs)
of Eq. (1) and T (xk ; xs + Dx, zs + Dz ) is expressed as

DT (xk ; xs , zs , Dx , Dz) = T (xk ; xs + Dx , zs + Dz) - T (xk ; xs , zs)
=

Dz
+
c0

(zs + Dz)2 + (xs + Dx - xk )2
c0

intervals of 44 ´ DTs¢ = 36.52 ms. Therefore, the instantaneous change in the 2D displacement in Fig. 18 cannot bemeasured using the conventional speckle tracking method.
However, a previous study had already attempted tomeasure the 2D displacement of the heart wall using parallel
receive beamforming to improve the temporal resolution to1020 Hz,15) which shows that the superiority of the proposed
method cannot be obtained with respect to time resolution. A
theoretical comparison of the performance of this method andthe proposed method is the subject of future work.

zs2 + (xs - xk )2
c0

.

(5 )

Using the following three assumptions
zs2  (xk - xs )2 ,(6 )∣ Dx ∣  ∣ xK / 2 ∣ ,(7 )∣ Dz ∣  zs ,(8 )the time difference distribution DT (xk ; xs, zs, Dx, Dz ) in
Eq. (5) can be approximated by the following linear function
regarding the element position xk as

SJ1040-10© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

(a)

(b)

Fig. 19. (Color online) (a) Amplitudes of the element signals with {T (xk ; xs , zs )} and {T (xk ; xs , zs ) + DT} acquired before and after displacing the positionof the probe by (Dx, Dz ) = (1 mm, 0.2 mm), respectively. (b) Magnitude distribution of the 2D complex cross-correlation ∣ Cn (Dx, Dz; xs , zs )∣.

Fig. 20. (Color online) Time difference distribution
DT (xk ; xs , zs, Dx, Dz ) for (Dx, Dz ) = (-1 mm, -0.2 mm). The blueline shows the true value obtained in Eq. (5), the green line shows theapproximated value in Eq. (9), and the red line shows the result by linearlyﬁtting Eq. (9) to the actual data of the blue line.

DT (xk ; x s, z s, Dx, Dz )»x2x1 ⎧⎛ x s + Dx
1 (x s + Dx )2+ s ⎞ xk +- s + 2 Dz ⎫.
⎬
c0 ⎨
zs ⎠
2 z s + Dz
2z s⎩⎝ z s + Dz
⎭
⎜

(9 )

⎟

(Dx, Dz ) = (-1 mm, -0.2 mm). The blue line represents
the true value obtained using Eq. (5), the green line
represents the approximated value in Eq. (9), and the redline shows the result obtained by linearly ﬁtting Eq. (9) to theactual data of the blue line. From the slope and intercept of
the red line, the 2D displacement estimated by the method of
xp , Dzp ) =a previous study28) was obtained as (D(-0.9 mm, - 0.14 mm). Compared to the previous study,
the 2D displacement estimated by the proposed method wascloser to the set value, as shown in Fig. 19.
In another previous study,29) after obtaining{DT (xk ; xs, zs, Dx, Dz )}, two elements were selected, andthe 2D displacement was estimated from the simultaneous
equations of the round-trip propagation paths. This method
yielded a lower signal-to-noise ratio than using all elements;
therefore, 2D displacement estimation is expected to be lessaccurate. Based on these comparisons, it can be concluded
that the proposed method can estimate the 2D displacement
more accurately than previous studies.28,29)
4. Conclusions

28)

In a previous study, the measured time difference distribution DT (xk ; xs, zs, Dx, Dz ) was linearly ﬁtted using the
xp , Dzp )least-squares method, and the 2D displacement (Dwas estimated from the slope and intercept of Eq. (9).
We compared the 2D displacements of the agar phantom
estimated by the proposed method and this previous study ina water tank experiment of Fig. 3 as follows. The phantom
was set in a water tank, and the probe was moved in the x–z
stage. The displacement of the probe was set as Δz = 0.2 mm
and Δx = 1 mm in the axial and lateral directions, respectively. The target was set in the phantom at position(xs, zs) = (-7 mm, 30 mm), assuming that the lateral position of the target xs moved away from the centre of the probe.
Figure 19(a) shows the amplitudes of the element signals
acquired before and after displacing the probe by(Dx, Dz ) = (1 mm, 0.2 mm), respectively. The time distributions {T (xk ; xs, zs)} and {T (xk ; xs, zs) + DT} are indicated by the red lines in Fig. 19(a). For the assumed position
x , zs + Dz )} and
x , Dz ), the distributions {T (xk ; xs + D(Dx , zs + Dz ) + DT} are indicated by the black{T (xk ; xs + D
lines in Fig. 19(a). Figure 19(b) shows ∣ Cn (Dx, Dz; xs, zs)∣
and the 2D displacement was estimated to bex , Dz ) = (-1.04 mm, -0.2 mm), which was consistent(Dwith the true value (-Dx, -Dz ) = (-1 mm, -0.2 mm).
Figure 20 shows the time difference distributionforDT (xk ; xs, zs, Dx, Dz )This paper proposed a novel 2D displacement estimation
method with a high temporal resolution using a 2D complex
cross-correlation comprising the element RF signals between
frames based on the received time distribution of scatteredwaves from the target. A phantom experiment was conducted
to examine the accuracy of the method, and the human heart
experiment veriﬁed its usefulness in the evaluation of cardiac
function.
In the phantom experiment, the estimated 2D displacementof the phantom was consistent with the set displacement,
conﬁrming that the proposed method can estimate the 2D
displacement of the phantom. In the in vivo experiment, wemeasured the 2D displacement of the IVS during one cardiac
cycle. However, during the rapid ﬁlling phase, the maximum
value of the cross-correlation decreased because of thechange in the RF waveforms. However, there are certain
aspects to be considered. The effect of changes in eachelement of the RF signal scattered from the myocardium
caused by contraction and dilation on the correlation should
be considered to further improve the accuracy. Moreover, in
this study, the length (ΔT) and width (K + 1) of thecorrelation window were ﬁxed. Thus, to improve theaccuracy, these parameters should be considered.
Furthermore, the limitations, that is, the lower and higher
limits, in estimating the instantaneous 2D displacement by
the proposed method should be considered.

SJ1040-11© 2023 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 62, SJ1040 (2023)

K. Hara et al.

Acknowledgments

This work was supported in part by JSPS KAKENHI
21H03835.

1) W. N. McDicken, G. R. Sutherland, C. M. Moran, and L. N. Gordon,
Ultrasound Med. Biol. 18, 651 (1992).
2) H. Kanai, H. Hasegawa, N. Chubachi, Y. Koiwa, and M. Tanaka, IEEE
Trans. Ultrason. Ferroelectr. Freq. Control 44, 752 (1997).
3) A. Heimdal, A. Støylen, H. Torp, and T. Skjærpe, J. Am. Soc. Echocardiogr.
11, 1013 (1998).
4) G. R. Sutherland, G. Di Salvo, P. Claus, J. D’hooge, and B. Bijnens, J. Am.
Soc. Echocardiogr. 17, 788 (2004).
5) J. Grondin, V. Sayseng, and E. E. Konofagou, IEEE Trans. Ultrason.
Ferroelectr. Freq. Control 64, 1212 (2017).
6) Y. Obara, S. Mori, M. Arakawa, and H. Kanai, Jpn. J. Appl. Phys. 60,
SDDE02 (2021).
7) P. Liu, J.-W. Muller, H.-M. Schwab, and R. Lopata, Front. Phys. 10, 1(2022).
8) D. Durrer, R. Th. Van Dam, G. E. Freud, M. J. Janse, F. L. Meijler, and R.
C. Arzbaecher, Circulation 41, 899 (1970).
9) H. Kanai, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52, 1931 (2005).
10) A. Hayashi, S. Mori, M. Arakawa, and H. Kanai, Jpn. J. Appl. Phys. 58,
SGGE05 (2019).
11) A. Basarab, P. Gueth, H. Liebgott, and P. Delachartre, IEEE Trans. Ultrason.
Ferroelectr. Freq. Control 56, 945 (2009).
12) Y. Honjo, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 49, 07HF14(2010).
13) D. Asari, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 53, 07KF21(2014).
14) M. V. Andersen, C. Moore, K. Arges, P. Søgaard, L. R. Østergaard, S.
E. Schmidt, J. Kisslo, and O. T. Von Ramm, Ultrasond Med. Biol. 42, 2606
(2016).

15) P. Joos, J. Porée, H. Liebgott, D. Vray, M. Baudet, J. Faurie, F. Tournoux,
G. Cloutier, B. Nicolas, and D. Garcia, IEEE Trans. Ultrason. Ferroelectr.
Freq. Control 65, 720 (2018).
16) M. Mozumi, R. Nagaoka, and H. Hasegawa, Jpn. J. Appl. Phys. 58,
SGGE02 (2019).
17) M. Orlowska, A. Ramalli, A. Petrescu, M. Cvijic, S. Bézy, P. Santos,
J. Pedrosa, J.-U. Voigt, and J. D’hooge, IEEE Trans. Ultrason. Ferroelectr.
Freq. Control 67, 1764 (2020).
18) H. Chen, T. Varghese, P. S. Rahko, and J. A. Zagzebski, Ultrasonics 49, 98(2009).
19) V. Mor-Avi et al., J. Am. Soc. Echocardiogr. 24, 277 (2011).
20) J.-U. Voigt et al., Eur. Heart. J. Cardiovasc. Imaging 16, 1 (2015).
21) D. P. Shattuck, M. D. Weinshenker, S. W. Smith, and O. T. von Ramm, J.
Acoust. Soc. Am. 75, 1273 (1984).
22) S. Langeland, J. D’hooge, H. Torp, B. Bijnens, and P. Suetens, UltrasoundMed. Biol. 29, 1177 (2003).
23) H. Hasegawa and H. Kanai, J. Med. Ultrason. 38, 129 (2011).
24) C. Papadacci, M. Pernot, M. Couade, M. Fink, and M. Tanter, IEEE Trans.
Ultrason. Ferroelectr. Freq. Control 61, 288 (2014).
25) R. Nagaoka, M. Mozumi, and H. Hasegawa, Jpn. J. Appl. Phys. 58,
SGGE10 (2019).
26) N. Furusawa, S. Mori, M. Arakawa, and H. Kanai, Jpn. J. Appl. Phys. 58,
SGGA08 (2019).
27) K. Sugahara, S. Mori, M. Arakawa, and H. Kanai, Jpn. J. Appl. Phys. 60,
SDDE05 (2021).
28) S.-L. Wang, M.-L. Li, and P.-C. Li, IEEE Trans. Ultrason. Ferroelectr. Freq.
Control 54, 70 (2007).
29) P. Kruizinga, F. Mastik, J. G. Bosch, N. de Jong, A. F. W. van der Steen,
and G. van Soest, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 62, 1733
(2015).
30) K. Hara, S. Mori, M. Arakawa, and H. Kanai, Proc. Symp. Ultrason. Electr.
43, 3Pa5-11 (2022).
31) W. Kroebel and K.-H. Mahrt, Acta Acustica United Acustica 35, 154 (1976).
32) C. Ma and T. Varghese, Ultrasonics 53, 782 (2013).
33) T. Asai, H. Taki, and H. Kanai, IEEJ Trans. Electron. Inf. Syst. 137, 573(2017).

SJ1040-12© 2023 The Japan Society of Applied Physics