Japanese Journal of Applied Physics 49 (2010) 07HF14REGULAR PAPER

Two-Dimensional Tracking of Heart Wall for Detailed Analysis
of Heart Function at High Temporal and Spatial Resolutions
Yasunori Honjo1 , Hideyuki Hasegawa1;2 , and Hiroshi Kanai2;1
12Graduate School of Biomedical Engineering, Tohoku University, Sendai 980-8579, Japan
Graduate School of Engineering, Tohoku University, Sendai 980-8579, Japan

Received November 24, 2009; accepted March 15, 2010; published online July 20, 2010
For noninvasive and quantitative measurements of global two-dimensional (2D) heart wall motion, speckle tracking methods have been developed
and applied. These methods overcome the limitation of tissue Doppler imaging (TDI), which is susceptible to aliasing, by directly tracking
backscattered echoes by pattern matching techniques, i.e., the cross-correlation or the sum of absolute differences, in real time. In theseconventional methods, the frame rate (FR) is limited to about 200 Hz, corresponding to the sampling period of 5 ms. However, myocardial function
during the isovolumic contraction period obtained by these conventional speckle tracking methods remains unclear owing to low temporal andspatial resolutions of these methods. Moreover, the accuracy of the speckle tracking method depends on an important parameter, i.e., the size ofthe correlation kernel. To track backscattered echoes accurately, it is necessary to determine the optimal kernel size. However, the optimal kernel
size has not been thoroughly investigated. In this study, correlation kernel size, which determines the tracking accurately, was optimized byevaluating root mean squared (RMS) errors in the lateral and axial displacements of a phantom estimated by speckle tracking methods at highspatial and temporal resolutions. For this purpose, the RF data from the longitudinal-axis cross-sectional view for the interventricular septum (IVS)
were acquired on the basis of parallel beam forming (PBF) to improve temporal and spatial resolutions. A wide transmit beam scanned in 7 different
directions sparsely and 16 receiving beams were generated for each transmission. The RF data of the phantom and the heart wall were obtained athigh spatial (angle intervals of scan lines: 0.375 degrees) and temporal [frame rate (FR): 1020 Hz] resolutions. The determined optimal size of thecorrelation kernel was 7:9  4:8 mm. Estimated displacements of the phantom were in good agreement with the actual displacement at an RMS
error of 0.34 mm. Furthermore, the IVS motion during the isovolumic contraction (IC) was analyzed in detail. The speckle tracking method using the
optimal kernel size 7:9  4:8 mm was applied to multiple points in IVS to estimate 2D displacements during the IC period. In this period, a rapid
displacement of IVS at a small amplitude of 1.5 mm, which suggests the expansion of the left ventricle and has not been measured by conventional
tracking methods at a low temporal resolution, was detected by 2D tracking. Furthermore, the displacement on the apical side was found to bedelayed by 10 ms compared with that on the basal side. These results indicate the potential of this method in the high-accuracy estimation of 2D
displacements and detailed analyses of physiological function of the myocardium. # 2010 The Japan Society of Applied Physics
DOI: 10.1143/JJAP.49.07HF14

1.

Introduction

For noninvasive and quantitative measurement of global
two-dimensional (2D) displacements, speckle tracking methods have been developed and applied.1–3) These methods
overcome the limitation of tissue Doppler imaging (TDI),4–8)
which is susceptible to aliasing, by directly tracking backscattered echoes using regional pattern matching techniques,i.e., the cross-correlation function or the sum of the absolute
diﬀerences, in real time.9,10) However, the frame rate (FR) in
conventional speckle tracking methods is limited to about
200 Hz, corresponding to the sampling period of heart wallmotion of 5 ms. In the isovolumic contraction (IC) period,
the heart wall moves rapidly during a short period of about
10 ms.11) Therefore, the detailed analysis of the heart
function has been limited owing to a low temporal resolutionof conventional speckle tracking methods. The continuous
observation of the heart wall motion requires a temporal
resolution much shorter than 10 ms.
Wang et al. proposed a method for multi-sector ultrasound imaging with retrospective electrocardiogram (ECG)
triggering, which achieved a high frame rate of 481 Hz.12)
Alternatively, the required high temporal resolution (FR >500 Hz) was realized by scanning the heart wall sparsely atthe expense of lateral spatial resolution.13–15) Therefore,
the spatial resolution in the measurement of displacementin the lateral direction has been signiﬁcantly degraded andthe speckle tracking method cannot be used with sparse
scanning. In this study, the ultrasonic RF data were acquiredon the basis of parallel beam forming so that the 2D
displacements of the heart wall and the myocardial functionE-mail address: honjo@us.ecei.tohoku.ac.jp

are measured in detail at a high temporal resolution less than
1 ms (corresponding to a frame rate of 1020 Hz) without
signiﬁcant degradation of spatial resolution.
Moreover, the accuracy of the speckle tracking method
depends on an important parameter, i.e., the size of thecorrelation kernel. To track backscattered echoes accurately,
it is necessary to determine the optimal kernel size.
However, the optimal kernel size has not been thoroughly
investigated. Trahey et al. reported that blood ﬂow can bemeasured by tracking echoes from the popliteal vein bynormalized correlation with a correlation kernel size of
2:0  1:5 mm2 .16) Similarly, Bohs et al. accurately determined the velocity of a string phantom, which moves twodimensionally at velocities up to 2.5 m/s, with a correlation
kernel size of 1:0  2:0 mm2 . The blood ﬂow in the lateral
direction was qualitatively imaged in in vitro experiments
using a correlation kernel size of 0:8  0:8 mm2 .17) However, in these studies, correlation kernel size has not been
optimized. Therefore, in this study, the optimal kernel size
was determined in a phantom experiment and then, thedetermined optimal kernel size was applied to the in vivo
measurement of 2D displacements of the heart wall by
normalized cross-correlation between RF echoes.
2.

Principles

2.1 Acquisition of RF signals at high frame rate

In conventional sector scanning, scan lines are generated
by focusing ultrasound to a certain spatial point in bothtransmission and reception with properly selected delays forthe signals applied to and received by transducer elements.
As illustrated in Figs. 1(a-1) and 1(a-2), the directions of
the transmit and receive beams coincide. In Figs. 1(a-1) and1(b-1), transmit beams are shown by dashed lines and in

07HF14-1# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14(a–1)

Y. Honjo et al.

(a–2)

(a–1)

(b–1)
top

sector probe

10 mm

sector probe

...

...

10 mm

15 mm
bottom

δθ

...

...

...

...

N tr =7

1
(a–2)

δθ

apical
side

(b–1)

(b–2)

...
sector probe

...

...

RV

...

apical
side

IVS
LV

basal
side

(c) RF acquisition area

10 mm

(RV)
right ventricular

RV
IVS
LV

basal
side

IVS
left ventricular(LV)

sector probe...

N t =7

(b–2)

10 mm

...

Fig. 2. Comparison of cross-sectional images obtained by sparse
scanning (a) with parallel beam forming (b). These are cross-sectional
images of a phantom (1) and interventricular septum (IVS) (2). (c)Acquisition area of RF signals in the longitudinal-axis view.

1
N r . δθ...
...

δθ =0.375

N r =16

Fig. 1. Transmission and reception of ultrasonic beams. (a) Conventional sector scanning. (b) Parallel beam forming. (1) Transmission.
(2) Reception.

Figs. 1(a-2) and 1(b-2), receive beams are shown by solidlines, respectively. FR in conventional scanning is given by
dividing the pulse repetition frequency fPRF by the number
Ntr of transmit beams, i.e., FR ¼ fPRF =Ntr . As shown in
Figs. 2(a-1) and 2(a-2), on the basis of the conventional
sector scan, Ntr ¼ 7 ultrasonic beams scanned a phantom anda heart wall sparsely at angular intervals of  ¼ 5:6 and ahigh temporal resolution of about 1/600 s using a 3.75 MHz
sector-type probe of ultrasonic diagnostic equipment (Aloka
SSD-6500). However, owing to the sparse scanning, thelateral spatial resolution is signiﬁcantly degraded.
In contrast, parallel beam forming (PBF)18,19) allows
receive beam forming in several directions for each ultrasound transmission. In this study, as shown in Figs. 2(b-1)
and 2(b-2), the ultrasonic beams scanned densely on thebasis of PBF so that the detailed analysis of the heart wallmotion is realized in a wide area at a high temporal
resolution of about 1 ms without signiﬁcant degradation of
the lateral spatial resolution. As shown in Fig. 1(b-1), in
this study, Nt ¼ 7 plane waves were transmitted in sevendiﬀerent directions at angular intervals of 6 at fPRF of
7170 Hz using a 3.75 MHz sector-type probe. In Fig. 1(b-2),
for each transmission, Nr ¼ 16 receiving beams weregenerated at angular intervals of  ¼ 0:375 . The frame
rate obtained by PBF is described by FR ¼ fPRF =Nt ¼1020 Hz. In addition, the total number of received beams
was Nt  Nr ¼ 112 and their interval (the lateral angle
spacing) was  ¼ 0:375 . The sampling frequency fs of the
RF signal was 15 MHz. In this study, as described above, ahigh FR (1020 Hz) acquisition of RF echoes was achieved ata spatial resolution much higher than that of the sparse scanwith a wide lateral view of Nt  Nr   ¼ 42 .
2.2 Evaluation of the normalized cross-correlation
function between RF echoes

The algorithm based on the correlation function between

ultrasonic RF echoes is one of the most common methods
to estimate the heart wall motion.3,20–22) In this study, 2D
displacements of the heart wall are estimated by employing
the normalized cross-correlation function between RF
echoes. Figure 3 shows the geometry of speckle tracking
based on the regional pattern matching technique using thenormalized cross-correlation function (2D tracking). The
center of a correlation kernel in the n-th frame is assigned
at the lateral and axial coordinates ð; dÞ  ðm; kÞ [point]( ¼ m  ; d ¼ k  d). The size of the correlation kernelin the lateral direction, Wl , is deﬁned by multiplying
the number of beams, 2Nl þ 1 ðNl Þ [points] and thelateral angle spacing  (deg/point) [lateral spacing l(mm/point) of the sampled RF signals] as Wl ¼ ð2Nl þ 1Þ   ð2Nl þ 1Þ  l. In this study, the lateral angle spacing
 was set to 0.375 as described in §2.1. The lateralspacing l (mm/point) at a depth of d was approximately
given by l  2d  sin:
ð1Þ
2
The size of the correlation kernel in the axial direction, Wd ,is deﬁned by multiplying the number of sampled RF signals,
2Nd þ 1 ðNd Þ [points], by the axial spacing d (mm/point)
as Wd ¼ ð2Nd þ 1Þ  d (mm). The axial spacing d is
deﬁned by the sampling frequency fs and the sound speed c0of the medium as follows:
c0
d ¼
:
ð2Þ
2  fsSpeckle size depends on the width of the ultrasonicbeam and pulse length.23,24) Therefore, in this study, thesizes of the correlation kernel in the lateral direction Wl ¼
ð2Nl þ 1Þ  l [corresponding to angular width ð2Nl þ 1Þ  ]and that in the axial direction Wd ¼ ð2Nd þ 1Þ  d weredeﬁned using the widths l and d at half maxima of
the lateral proﬁle of an ultrasonic ﬁeld and the envelopeof an ultrasonic pulse at a depth of 50 mm (focal area),
respectively, as shown in Fig. 4. The shape of thecorrelation kernel was a bivariate normal probability
distribution. The variances l2 (mm) and d2 (mm) of thenormal distribution in the lateral and axial directions were
determined as

07HF14-2# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14

Y. Honjo et al.

m +2

n –th framem . δθ

m+1 beam m

m –1

m –2

σl
2 σl

dn

10 mm

2 σl
Nl

Nl

Wd

m –2

σd

m –1

2σ d

Nd

Wl

(m , k )correlation kernel(m , k )

2σ d

Nd
k

(m , k )m th beam

m +6

m+5

m +4

m +3

m +2

m +1 beam m( n +Δ NF )–th frame(m + Δ m n , k + Δ k n)
(m , k )

k–3

Δk n

k
(m + Δ m n , k + Δ k n)best matching( m , k)Δm nFig. 3.

(Color online) Calculation of normalized cross-correlation function between RF signals in the n-th and ðn þ NF Þ-th frames.

(a–2) 1
0.9

(a–1)

normalized amplitude

stage
controller
ultrasonic probe

water
50 mm
wire (35 μ m)

0.8
0.7
0.6
0.5
1.6 mm

0.4
0.3
0.2
0.1

water tank

automatic X stage

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6lateral distance [mm](b–2) 1
0.9

ultrasonic probe

normalized amplitude

(b–1)

0

water
50 mm
wire (35 μ m)

0.8
0.7
0.6
0.5
1.1 mm

0.4
0.3
0.2
0.1
0

water tank00.5 1 1.5 2 2.5 3 3.5 4 4.5
distance in the direction of depth [mm]5Fig. 4. (a) Schematic image of system for measurement of beam width and measured beam profile. (b) Schematic image of system for measurementof pulse length and measured pulse envelope. The depth d was set at 50 mm (focal area).

l ¼   l
d ¼   d

(mm);
(mm);

ð3Þ
ð4Þ

where  is a constant. The correlation kernel wði; jÞ applied
to the RF signal rf n ði; jÞ in the n-th frame at i-th lateral point
and j-th depth point was deﬁned as follows:
07HF14-31wði; jÞ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ exp
2 detðÞ iDði; jÞ ¼
;
j




1
T
1
Dði; jÞ  Dði; jÞ ; ð5Þ
2
ð6Þ# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14

¼

l2

0

0

d2

!
¼

Y. Honjo et al.

ð  lÞ200ð  dÞ2

!
;

ð7Þ

where wði; jÞ is a window function corresponding to thebivariate normal probability distribution, Dði; jÞ is a vector
from the center of the correlation kernel to the coordinate ði; jÞ (i ¼ Nl ; Nl þ 1; . . . ; 1; 0; 1; . . . ; Nl ; j ¼ Nd ;Nd þ 1; . . . ; 1; 0; 1; . . . ; Nd ), and  1 is the inverse
covariance matrix of the bivariate normal probability
distribution . In this study, the sizes Wl and Wd of thecorrelation kernel in the lateral and axial directions weredeﬁned by 2ðl ; d Þ ¼ 2ð  l;   dÞ so that RF
signals in the correlation kernel are included in the 95%
region of the bivariate normal probability distribution.
The position of the correlation kernel with a size ofðWl  Wd Þ ¼ ðð2Nl þ 1Þ    ð2Nd þ 1Þ  dÞ is assigned
manually in the ﬁrst frame and its central position is tracked
in the successive frames. The normalized cross-correlation
function r;D;n ðmn ; kn ; NF Þ at the lateral shift mn(points) and axial shift kn (points) is calculated from RF
signals frf n ðm; kÞg and frf nþNF ðm þ mn ; k þ kn Þg (m:beam number, k: sampled point number of RF signals in
the axial direction) in the n-th and ðn þ NF Þ-th frames asfollows:

¼

kþN
Xd

1
fwði
A i¼m Nl j¼k Nd

R;D;n ðm0n ; kn0 ; NF Þ
("
N
Nl =2
d =2
X
X
¼

r;D;n ðmn ; kn ; NF Þ

mn ¼Nl =2

kn ¼ Nd=2#)

 sincðl Þ  wh ðl Þ  sincðd Þ  wh ðd Þ ;
l ¼
d ¼

m0n L0l

mn LlLlkn0 L0d

kn Ld
Ld

 ;

ð11Þ

 ;

ð12Þ

sin 
;
ð13Þ

wh ðÞ ¼ 0:423 { 0:498  cosð2Þ þ 0:0792 cosð4Þ; ð14Þr;D;n ðmn ; kn ; NF ÞmþN
Xl

normalized cross-correlation function are not suﬃcient for
the estimation of 2D displacements.
To overcome this problem, the above normalized crosscorrelation function r;D;n ðmn ; kn ; NF Þ shown in eq. (8)was interpolated by reconstructive interpolation20) to improve the accuracy in estimation of the displacement. Thismethod is a space-domain technique based on the signal
reconstruction theory using the sinc function, that is, thesampled normalized cross-correlation function R;D;n ðm0n ;kn0 ; NF Þ at an arbitrary lateral shift m0n and axial shiftkn0 is given by the following low-pass ﬁltering after
increasing spatial sampling frequency:sincðÞ ¼
m; jkÞ  rf n ði; jÞ wði m; j kÞ  rf nþNF ði þ mn ; j þ kn Þg; ð8Þ
A ¼ ðm;kÞ  ðmþmn ;kþkn Þ ;ð9Þ
vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
u mþNl kþNdu X X
ðm;kÞ ¼ t
w2 ði m; j kÞ  rf 2n ði; jÞ; ð10Þi¼m Nl j¼k Nd

where ðm;kÞ is the standard deviation of the windowed RF
signals frf n ðm; kÞg in the correlation kernel with a size ofðWl  Wd Þ ¼ ðð2Nl þ 1Þ    ð2Nd þ 1Þ  dÞ and NF is
the frame interval for calculation of correlation. As
illustrated in Fig. 3, RF signals in the correlation kernel in
the n-th frame are compared with those in the ðn þ NF Þ-thframe. The lateral and axial displacements are determinedfrom the lateral shift c
mn and axial shift kbn , which
maximize the normalized cross-correlation function
fr;D;n ðmn ; kn ; NF Þg.
2.3 Reconstructive interpolation

The resolutions in the estimation of the lateral and axial
displacements, which are determined from the lateral shiftcmn and axial shift kbn , depend on the lateral angle spacing (deg) [lateral spacing l (mm)] and the axial spacing d(sampling interval of the RF signals) of the correlation
function. In this study,  and d were set at 0.375 and51.3 mm (sampling frequency fs ¼ 15 MHz), respectively.
The lateral spacing l at a typical depth of the interventricular septum (IVS) (d ¼ 50 mm) is 330 mm and it is muchlarger than the axial spacing d of 51.3 mm. In this study, the2D displacements of the heart wall between consecutive
frames are expected to be very small because it is measured
at a high temporal resolution T of 1 ms. Therefore, thesampling intervals  ð lÞ and d of the coarsely sampled

where Nl and Nd are the numbers of the original normalized
cross-correlation function in the lateral and axial directions
used for interpolation, respectively, Ll , Ld , and L0l , L0d are theintervals of normalized cross-correlation functions in thelateral and axial directions before and after interpolation,
respectively, sincðÞ is the sinc function, and wh ðÞ is theBlackman window.
The lateral and axial average velocities, denoted byvl ðn; ; dÞ and vd ðn; ; dÞ, between the n-th and ðn þ NF Þ-th
c0 n andframes are respectively given by the lateral shift m0caxial shift k n , determined so that the normalized crosscorrelation function in eq. (11) takes the maximum asfollows:
c0 n  L0l
mvl ðn; ; dÞ ¼
(mm/s);
ð15Þ
NF  T
0  L0kcnd(mm/s);
ð16Þ
vd ðn; ; dÞ ¼
NF  T
Then, the lateral displacement xl ðn þ 1; ; dÞ and the axial
displacement xd ðn þ 1; ; dÞ at the coordinate ð; dÞ in theðn þ 1Þ-th frame are estimated as follows:xl ðn þ 1; ; dÞ ¼ xl ðn; ; dÞ þ vl ðn; ; dÞ  T (mm); ð17Þ
xd ðn þ 1; ; dÞ ¼ xd ðn; ; dÞ þ vd ðn; ; dÞ  T (mm): ð18ÞThe instantaneous position of the correlation kernel in theðn þ 1Þ-th frame is determined by xl ðn; ; dÞ and xd ðn; ; dÞ.
3.

Determination of the Optimal Size of a Correlation
Kernel by Basic Experiment

3.1 Experimental system

Figure 5(a) shows an experimental system. The motion ofa phantom was measured to determine the optimum size
cd ) of the correlation kernel. As shown in Fig. 5(b),(cWl  W

07HF14-4# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14

stage
controller

Y. Honjo et al.

lateralaxial(a)d =50 mm(b)

ultrasonic probe

water
automatic X–Y stage

water tank

=
1

top

...

300
...

10 mm

silicone rubber

j

=

k

15 mm
500
Md
m =70 ... 40
bottom( i = M l ... 1)
area of 2D tracking

Fig. 5. (Color online) (a) Schematic image of system for basic experiment. (b) Acquired cross-sectional image of silicone rubber. The 2D
displacements of the silicone rubber in the area were estimated. The beam number m for the displacement estimation was set at intervals of five beams
from 40 to 70. The depth k was set at intervals of 50 points from 300 to 500 from the top of the silicone rubber.

the phantom was made from a silicone rubber (Momentive
Performance Materials TSE3503), containing 5% graphite
powder by weight to provide suﬃcient scattering. The
phantom was a plate (sound speed c0 : 1000 mm/s) with athickness of 15 mm. The motion velocity of the silicone
rubber was controlled by an automatic X–Y stage. The
phantom was moved two-dimensionally. In this study, since
the instantaneous rotation component of the regional heart
wall motion is negligibly small at a high temporal resolution,
this experiment excluded rotational motion. Both the lateralvelocity vl0 ðn; ; dÞ and axial velocity vd0 ðn; ; dÞ were set
at a constant speed of 5 mm/s. However, the velocities
vl0 ðn; ; dÞ and vd0 ðn; ; dÞ of the silicone rubber in the lateraland axial directions were much lower than the maximum
velocities of the actual heart wall motion (about 70 mm/s
for both) obtained in separate in vivo experiments by the1D phased-tracking method.11,25,26) Therefore, 2D displacements of silicone rubber between two frames (NF ) in thelateral and axial directions, xl ðn; ; dÞ and xd ðn; ; dÞ, wererespectively set to be similar to those of the heart wall by
increasing the frame interval NF for calculating correlation
function as follows. The frame interval NF was set at 25(corresponding 24 ms) in the basic experiment so that the
2D displacements of the silicone rubber between frames(NF  T) correspond to the fast motion of the heart wall
between two consecutive frames, corresponding to 0.98 ms(NF ¼ 1) in in vivo experiments.
3.2 Determination of the optimal size of the correlation
kernelEstimated displacement depends on the size of the correlation kernel. Therefore, determination of the optimal correlation kernel size is important. When the size of thecorrelation kernel is set to be small, there may be several
regions in the post frame, which have echo patterns similar
to those in the correlation kernel in the previous frame.
These multiple regions, which have similar echo patterns,
introduces misestimation of displacement.22) On the other
hand, when the size of the kernel is set to be large, theecho pattern in the correlation kernel is more unique andthe displacement estimation is less susceptible to noise.
However, a proﬁle around a peak of the normalized crosscorrelation function becomes wider because it is rare thatthe echo patterns in two frames coincide, which degrades
the accuracy in the estimation and the spatial resolution.17)
In this study, the sizes Wl and Wd of the correlation kernelin the lateral and axial directions were deﬁned by 2lð¼ 2  lÞ and 2d ð¼ 2  dÞ, respectively. To
cl  Wcd ) of the correlation
determine the optimal size (Wkernel, the 2D displacements estimated using diﬀerent
values of  were examined using the root-mean-squared(RMS) error of the estimated displacement from the actual
displacement. The RMS error in the estimated 2D displacements, ðÞ, which was evaluated for each value of , wascalculated from the lateral displacement xl ð; d; Þ and axial
displacement xd ð; d; Þ estimated for each position ð; dÞi ðm  ; k  dÞi [as shown in Fig. 5(b), i ¼ 1; 2; . . . ; Ml  Md
(m ¼ 40; 45; . . . ; 65; 70-th beam, k ¼ 300; 350; . . . ; 450; 500th sampling point)] and the actual lateral displacement xl0 ðnÞ
and axial displacement xd0 ðnÞ as follows:
vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃuMNl Md Xu 1
XðÞ ¼ t½xl ðn; i; Þ2 þ xd ðn; i; Þ2 ;Ml Md N i¼1 n¼1
ð19Þ
xl ðn; i; Þ ¼ xl ðn; i; Þ

xl0 ðnÞ;ð20Þ

xd ðn; i; Þ ¼ xd ðn; i; Þ

xd0 ðnÞ;

ð21Þ

where Ml and Md are the numbers of measured points inthe lateral and axial directions, respectively, and N is thenumber of averaged frames. In this experiment, as illustrated
in Fig. 5(b), the displacements at the positions M ð¼ Ml 
Md ¼ 7  5Þ ¼ 35 in the area of 2D tracking were measured
by normalized cross-correlation functions using diﬀerent
values  ¼ 0:2 { 2:0, corresponding to the kernel sizes Wl ¼
1:2 {12:9 mm and Wd ¼ 0:9 { 9:0 mm at a depth of 50 mm.
The original sampled normalized cross-correlation functionr;D;n ðm; k; NF Þ in eq. (8) was interpolated by factors
L0l =Ll ¼ 60 in the lateral direction and L0d =Ld ¼ 10 in theaxial direction, to obtain the cross-correlation function
R;D;n ðm0n ; kn0 ; NF Þ in eq. (11) with similar resolutionsin the lateral and axial directions.
Figure 6 shows the distribution of RMS errors, ðÞ,estimated using diﬀerent values of  of the correlation
kernels. When 2D displacements were estimated using asmall kernel size, the RMS error was large, as shown in

07HF14-5# 2010 The Japan Society of Applied Physics

Y. Honjo et al.

1.8

(a)

1.2

2 σl

2 σl
Wd

2 σd

σl

1.4

2 σd

1

Wl

0.8
0.6
0.4
0.2
0.2

0.4

1.2

0.9

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.5

3.8
5.1
6.4
7.7
9.0
10.3 11.6
kernel size in the lateral direction ( Wl ) [mm]

12.9

1.8

2.7
3.6
4.5
5.4
6.3
7.2
kernel size in the depth direction ( Wd ) [mm]9.0displacement xd ( n; θ, d ) in the axial direction [mm]this ﬁgure. This result shows that it is impossible to trackthe 2D displacements accurately using a small and overlarge
kernel sizes, such as ðWl  Wd Þ ¼ ð1:2  0:9Þ mm ( ¼ 0:2)
and ðWl  Wd Þ ¼ ð12:9  9:0Þ mm ( ¼ 2:0). On the other
hand, when the kernel size was set at ðWl  Wd Þ ¼ ð7:0 
4:8Þ mm ( ¼ 1:1), the RMS error was minimized tocl  Wcd ) of the correlation
0.34 mm. This optimal size (Wkernel was four times the widths of the ultrasonic beam andthe pulse length at a depth of 50 mm (focal area). In thisstudy, the lateral width Wl of sampled points was assumed
to be constant because the change in the lateral spacingl of sampled points in the region of interest (IVS, typical
depth and thickness: about 50 and 10 mm) is negligible
compared with the ultrasonic beam width l. The size ofcl ¼
the optimal correlation kernel in the lateral direction, W
7:0 mm at the depth d ¼ 50 mm, corresponds to 7.9 [thenumber of the beam in the kernel is 2Nl þ 1 ¼ 21, as shownin eq. (1)].

(1.0 x 2.0) mm(0.8 x 0.8) mm
actual

3.5
3
2.5
2
1.5
1
0.5
0

(b)

Fig. 6. Distribution of RMS errors ðÞ for difference sizes of kernels
ðWl  Wd Þ ¼ ½ð1:2 {12:9Þ  ð0:9 { 9:0Þ mm.

(7.0 x 4.8) mm(2.0 x 1.5) mm

4

–0.5

α

8.1

4.5

displacement xl ( n; θ, d ) in the lateral direction [mm]

1.6

σd

RMS error ε(α) in displacement estimation [mm]Jpn. J. Appl. Phys. 49 (2010) 07HF14

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6time [s]4.5(7.0 x 4.8) mm(2.0 x 1.5) mm

4(1.0 x 2.0) mm(0.8 x 0.8) mm
actual

3.5
3
2.5
2
1.5
1
0.5
0
–0.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6time [s]Fig. 7. (Color online) (a) Estimated displacement xl ðn; ; d Þ and actual
displacement xl0 ðnÞ in the lateral direction. (b) Estimated displacementxd ðn; ; d Þ and actual displacement xd 0 ðnÞ in the axial direction. These
displacements were estimated using various kernel sizes (Wl  Wd ) of
7:0  4:8, 1:0  2:0,16) 2:0  1:5, and 0:8  0:8 mm2 .17) The depth d ofeach kernel was set as 50 mm at the initial time.

(a) RF acquisition area

d =50 mm
RV
IVS

10 mm
Ao

LV
LA

3.3 Estimation of 2D motion of the silicone rubber using
the optimal kernel size

The 2D displacements of the silicone rubber were estimatedcd )by the 2D tracking using the optimal kernel size (cWl  Wof ð7:9  4:8 mmÞ  ð7:0  4:8Þ mm at the depth of d ¼
50 mm. Figures 7(a) and 7(b) show the lateral and axial
displacements obtained by the 2D tracking. The displacements estimated for M ¼ 35 points in the area of 2D
tracking are shown by red lines. The displacements in thelateral and axial directions with the correlation kernel size of(7:9  4:8 mm) were in good agreement with the actual
displacements (black line). On the other hand, the displacements estimated using the previously examined kernel
sizes,16,17) 2:0  1:5, 1:0  2:0, and 0:8  0:8 mm2 at thedepth d ¼ 50 mm in the initial frame, were markedly
diﬀerent from the actual displacement determined by thecontrol of the stepping motor. This result shows that the 2D
displacements were estimated accurately by the normalized
cross-correlation function with the determined optimal
cl  W
cd ).
correlation kernel size (W

(b)

PW

apical IVS
side
beam
80th

RV
...

basal
side

30th

Fig. 8. (Color online) (a) Acquisition area of RF signals from thelongitudinal-axis view. (b) The points employed for estimating 2D
displacements are shown by the red points on the cross-sectional imageof IVS. The beam numbers were set at 30 – 80 and the depth was set at
6.0 mm from the chest.

4.

In vivo Experimental Results

4.1 Acquisition of RF signals in IVS

Figure 8 shows a typical cross-sectional image (transthoracic parasternal longitudinal-axis view) of the heartobtained from a healthy 22-year-old male. Figure 8(a)shows the RF acquisition area. In the acquisition of RF
signals, the ultrasonic beam scanned Nt  Nr ¼ 112 diﬀerent directions densely to maintain a high spatial resolution. The frame rate was set at 1020 Hz. As shown in

07HF14-6# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14
apical area beam number basal area
80 75 70 65 60 55 50 45 40 35 30

10

ECG

PCG

0

20

time [ms]

30
40
50
60

80
90
100

apical side28

basal side

24 20 16 12
8
40
distance from 30th beam [mm]1st heart sound70

mitral valve closureR–waveisovolumic contraction

mitral valve closureR–waveisovolumic contraction(a)Y. Honjo et al.

ECG
PCG(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)

0

10

20

30
5 mm

40 50 60
time [ms]1st heart sound70 80 90 100

scale of the displacement

(1)
10 mm

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

–4

mitral valve closureR–waveisovolumic contraction

(b)
ECG

1st heart sound–4
RV side
0

30
35
40
45
50
55
60
65
70
75
80

4
8
12
16
20
24
LV side0

distance from 30th beam [mm]beam number

PCG

28

10 20 30 40 50 60 70 80 90 100
time [ms]Fig. 9. Displacements of IVS on the left ventricular side in the lateral
direction (a) and in the axial direction (b) from R-wave to isovolumic
contraction period.

Fig. 8(b), the motion of a region in IVS of the left
ventricular (LV) side was estimated by 2D tracking. Somefeatures of speckle patterns, such as the size of a speckle,
depend on the point spread function produced by theultrasonic ﬁeld.23,24) In this study, as shown in Figs. 2(b-1)
and 2(b-2), the RF data were acquired using a 3.75 MHz
sector-type probe of ultrasonic diagnostic equipment,
which was the same as that used in the basic experiment
so that the features of speckle size was not changed.
cl  W
cd Þ ¼
Therefore, the size of a correlation kernel ðWð7:9  4:8 mmÞ  ð7:0  4:8Þ mm determined by the basic
experiment was also used in the in vivo experiment.
4.2 Estimating motions of the interventricular septum fromthe time of the R-wave to isovolumic contraction period

Figure 9 shows the displacements in the lateral direction
xl ðn; ; dÞ and axial direction xd ðn; ; dÞ of IVS on the left
ventricular (LV) side plotted as functions of time fromthe R-wave during isovolumic contraction (IC) period.
Figures 10, 11, and 12 show the cross-sectional images of
IVS obtained at 3 ms intervals. The red lines in these images
show the direction and magnitude of motion of IVS.
As shown in Figs. 9(a) and 10, IVS on the LV side moved
to the apical side after the time of the R-wave. This motionof IVS in the apical area was larger than that of IVS in
the basal area approximately 20 ms after the R-wave. Thissuggests that IVS in the basal area was elongated. In

Fig. 10. (Color online) Motion of IVS after R-wave. The cross-sectional
images of IVS at every 3 ms (1)–(12). The red lines in these images show
the direction and magnitude of instantaneous displacement of IVS. The
yellow points in these images show the positions ð; dÞi of the i-th point ofinterest, which was set at intervals of five beams from 30 to 80.

addition, IVS in the basal area moved to the RV side after
the time of the R-wave. In contrast, IVS in the apical area
moved to the LV side. This is considered to result from anincrease in LV pressure owing to contraction of the apicalarea for the succeeding ejection phase. Figures 9(a) and 10
show that IVS moved from the apical side to the basal side
after 10 – 20 ms from the R-wave. These motions of IVS
propagated from the basal area to the apical area, as shownin Figs. 10(5)–(8). At the beginning of the isovolumic phase
after the mitral valve closed, IVS in the basal area moved
rapidly to the basal side. Moreover, IVS moved from the RVside to the LV side approximately 30 ms from the R-wave.
These motions of IVS propagated from the basal area to theapical area, as shown in Fig. 9(b). This suggests that IVS in
the basal area was pulled by the motion of the closure of themitral valve.
As shown in Figs. 9(b) and 11(1)–(5), IVS moved fromthe RV side to the LV side after 35 – 45 ms from the R-wave.
It is recognized that there is a time delay between the peak ofthe displacements in the basal area and that in the apicalarea. The displacement in the apical area was found to bedelayed by 10 ms from that in the basal area. Figures 9(b),11, and 12 show IVS further moved to the LV side after
45 ms. As shown in Figs. 9(a) and 12, IVS in the basal area

07HF14-7# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14

Y. Honjo et al.

mitral valve closureR–waveisovolumic contraction

mitral valve closureR–waveisovolumic contraction

ECG

ECG

PCG

PCG
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)

0

10

20

30
5 mm

40 50 60
time [ms]7080(1) (2) (3) (4) (5) (6) (7)90 1000scale of the displacement

10

20

30
5 mm

40 50 60
time [ms]7080

scale of the displacement

(1)
10 mm

(2)

(3)

(1)
10 mm

(2)

(3)

(4)

(5)

(6)

(4)

(5)

(6)

(7)

(8)

(9)

(7)

(10)

(11)

(12)

Fig. 12. (Color online) Motion of IVS from 70 to 90 ms after R-wave.
Cross-sectional images of IVS at every 3 ms (1)–(7). The red lines inthese images show the direction and magnitude of instantaneous
displacement of IVS. The yellow points in these images show thepositions ð; dÞi of the i-th point of interest, which was set at intervals of
five beams from 30 to 80.

Fig. 11. (Color online) Motion of IVS from 36 to 69 ms after R-wave.
Cross-sectional images of IVS at every 3 ms (1)–(12). The red lines inthese images show the direction and magnitude of instantaneous
displacement of IVS. The yellow points in these images show thepositions ð; dÞi of the i-th point of interest, which was set at intervals of
five beams from 30 to 80.

moved to the apical side at the time of 60 – 80 ms. Thesemotions of the IVS propagated from the basal area to theapical area.
It was found that IVS moved from the apical side to thebasal side and RV side during the isovolumic period toincrease LV pressure for the succeeding ejection phase. In
this study, the displacements in IVS during the IC period
could be observed at a high temporal resolution by the
proposed method.
5.

90 100similar to a typical heart wall motion. Furthermore, the IVS
motions during isovolumic contraction were analyzed indetail. In this period, a rapid displacement of IVS at a small
amplitude of 2.0 mm, which suggested the expansion of theleft ventricle and has not been measured by conventional
tracking methods at a low temporal resolution, was detectedby 2D tracking. Furthermore, the displacement on the apicalside was found to be delayed by 10 ms from that on the basalside. Although further investigation to conﬁrm whether theidentiﬁed relationship between the optimum kernel size andthe ultrasonic ﬁeld can be applied to other depths (6¼ 50 mm,
typical depth of IVS),27,28) these results indicate the potentialof this method in the high-accuracy estimation of 2D
displacements and the detailed analyses of the physiological
function of the myocardium.

Conclusions

In this study, an important parameter in 2D displacements
estimation was determined to be optimized by RMS errorsof the 2D displacements obtained using diﬀerent sizes ofcorrelation kernels at a high frame rate based on parallel
beam forming. The determined optimal size of the correlation kernel was 7:9  4:8 mm. The estimated displacementsof the silicone rubber were in good agreement with theactual displacement. The regional 2D displacements of theheart wall were measured at an extremely high temporal
resolution (frame rate about 1020 Hz) without signiﬁcant
increase in the lateral spacing of scan lines (beam intervals ¼ 0:375 ). The estimated displacements of IVS were

1) L. N. Bohs, B. J. Geiman, M. E. Anderson, S. C. Gebhart, and G. E.
Trahey: Ultrasonics 38 (2000) 369.
2) S. Langeland, J. D’hooge, H. Torp, B. Bijnens, and P. Suetens: UltrasoundMed. Biol. 29 (2003) 1177.
3) S. Langeland, J. D’hooge, T. Claessens, P. Claus, P. Verdonck, P.
Suetens, G. R. Sutherland, and B. Bijnens: IEEE Trans. Ultrason.
Ferroelectr. Freq. Control 51 (2004) 1537.
4) G. R. Sutherland, M. J. Stewart, K. W. Groundstroem, C. M. Moran, A.
Fleming, F. J. Guell-Peris, R. A. Riemersma, L. N. Fenn, K. A. Fox, andW. N. McDicken: J. Am. Soc. Echocardiogr. 7 (1994) 441.
5) P. Palka, A. Lange, A. D. Fleming, G. R. Sutherland, L. N. Feen, andW. N. McDicken: J. Am. Soc. Echocardiogr. 8 (1995) 659.
6) K. Miyatake, M. Yamagishi, N. Tanaka, M. Uematsu, N. Yamazaki, Y.
Mine, A. Sano, and M. Hirama: J. Am. Coll. Cardiol. 25 (1995) 717.

07HF14-8# 2010 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 49 (2010) 07HF14

Y. Honjo et al.

7) J. Gorcsan III, V. K. Gulati, W. A. Mandarino, and W. E. Katz: Am. HeartJ. 131 (1996) 1203.
8) A. Heimdal, A. Støylen, H. Torp, and T. Skjærpe: J. Am. Soc.
Echocardiogr. 11 (1998) 1013.
9) T. Numata, H. Hasegawa, H. Kanai, and M. Tanaka: Jpn. J. Appl. Phys.
46 (2007) 4900.
10) H. Yoshikawa, T. Azuma, K. Sasaki, K. Kawabata, and S. Umemura: Jpn.
J. Appl. Phys. 46 (2007) 4834.
11) H. Yoshiara, H. Hasegawa, H. Kanai, and M. Tanaka: Jpn. J. Appl. Phys.
46 (2007) 4889.
12) S. Wang, W. Lee, J. Provost, J. Luo, and E. Konofagou: IEEE Trans.
Ultrason. Ferroelectr. Freq. Control 55 (2008) 2221.
13) H. Kanai and Y. Koiwa: Ultrasound Med. Biol. 27 (2001) 481.
14) H. Kanai: IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 (2005)
1931.
15) H. Kanai: Ultrasound Med. Biol. 35 (2009) 382.
16) G. E. Trahey, J. W. Allison, and O. T. von Ramm: IEEE Trans. Biomed.
Eng. 34 (1987) 965.
17) L. N. Bohs, B. H. Friemel, B. A. McDermott, and G. E. Trahey:
Ultrasound Med. Biol. 19 (1993) 751.
18) D. P. Shattuck, M. D. Weinshenker, S. W. Smith, and O. T. von Ramm:

J. Acoust. Soc. Am. 75 (1984) 1273.
19) H. Hasegawa and H. Kanai: IEEE Trans. Ultrason. Ferroelectr. Freq.
Control 55 (2008) 2626.
20) I. Cespedes, Y. Huang, J. Ophir, and S. Spratt: Ultrason. Imaging 17(1995) 142.
21) F. Viola and W. F. Walker: IEEE Trans. Ultrason. Ferroelectr. Freq.
Control 50 (2003) 392.
22) D. Rappaport, D. Adam, P. Lysyansky, and S. Riesner: Ultrasound Med.
Biol. 32 (2006) 1181.
23) I. Akiyama, A. Hayama, M. Nakajima, S. Yuta, and T. Itoh: Denshi Joho
Tsushin Gakkai Ronbunshi D 71 (1988) 733 [in Japanese].
24) R. F. Wagner, S. W. Smith, J. M. Sandrik, and H. Lopez: IEEE Trans.
Sonics Ultrason. 30 (1983) 156.
25) H. Kanai, Y. Koiwa, Y. Saito, I. Susukida, and M. Tanaka: Jpn. J. Appl.
Phys. 38 (1999) 3403.
26) T. Kinugawa, H. Hasegawa, and H. Kanai: Jpn. J. Appl. Phys. 47 (2008)
4155.
27) G. E. Trahey and S. W. Smith: Ultrason. Imaging 10 (1988) 12.
28) S. W. Smith, G. E. Trahey, S. M. Hubbard, and R. F. Wagner: Ultrason.
Imaging 10 (1988) 29.

07HF14-9# 2010 The Japan Society of Applied Physics