Home

SearchCollections

Journals

About

Contact us

My IOPscience

Automated detection of arterial wall boundaries based on correlation between adjacentreceive scan lines for elasticity imaging

This content has been downloaded from IOPscience. Please scroll down to see the full text.
2015 Jpn. J. Appl. Phys. 54 07HF18(http://iopscience.iop.org/1347-4065/54/7S1/07HF18)
View the table of contents for this issue, or go to the journal homepage for more

Download details:
IP Address: 130.34.32.223
This content was downloaded on 08/08/2015 at 01:16

Please note that terms and conditions apply.

REGULAR PAPER

Japanese Journal of Applied Physics 54, 07HF18 (2015)
http://dx.doi.org/10.7567/JJAP.54.07HF18

Automated detection of arterial wall boundaries based on correlation
between adjacent receive scan lines for elasticity imaging
Yukiya Miyachi1,2, Hideyuki Hasegawa2,3, and Hiroshi Kanai3,2*
1

Medical Systems Research and Development Center R&D Management Headquarters, FUJIFILM Corporation, Kaisei, Kanagawa 258-8538, Japan
Graduate School of Biomedical Engineering, Tohoku University, Sendai 980-8579, Japan
3
Graduate School of Engineering, Tohoku University, Sendai 980-8579, Japan
E-mail: yukiya.miyachi@fujiﬁlm.com
2

Received December 12, 2014; revised March 12, 2015; accepted April 2, 2015; published online June 25, 2015
In our series of studies on the ultrasonic assessment of the regional elasticity of the arterial wall, the change in the thickness of the arterial wallcaused by the heartbeat was measured by the ultrasonic phased-tracking method. In the measurement of elasticity, the boundaries of the anterior
and posterior walls must be assigned in advance. Currently the boundaries are manually determined by the operator, which is time-consuming andresults in observer variability. In this paper, we propose an automated method for the detection of arterial wall boundaries using multiscale dynamic
programming, in which the cost function includes the correlation term between ultrasonic echoes in adjacent receive scan lines. The correlation
term enables boundary detection that is more robust than conventional methods against noises. The effectiveness of the proposed method wasvalidated using a phantom and applied to the in vivo measurements of carotid arteries. The root mean square error between the results obtainedby the proposed method and the manual assignment of the boundaries was signiﬁcantly improved. © 2015 The Japan Society of Applied Physics

1.

Introduction

Along with the recent changes in lifestyle in Japan, theincidence of coronary artery disease has increased.1) Atherosclerosis has been considered to be the main cause ofcoronary artery disease and cerebral infraction. Because there
are signiﬁcant diﬀerences in elasticity between normal andatherosclerotic arterial walls,2,3) the evaluation of arterial
wall elasticity is useful.4) Some markers that can be
noninvasively measured for the detection of atherosclerosis
have been identiﬁed, such as intima–media thickness (IMT)of the arterial wall, pulse wave velocity (PWV), and elasticityof the arterial wall.
Noninvasive ultrasonic B-mode imaging, where theamplitudes of received ultrasonic echoes are converted tobrightness, is an important method for studying the
progression and regression of atherosclerotic lesions in the
carotid artery.5–7) IMT is measured accurately from a B-mode
ultrasonic image and studies for more accurate measurementof IMT have been conducted.8–12) This measurement is used
as a sensitive indicator of atherosclerosis in population-based
studies. However, a conventional B-mode ultrasonic image
provides only information on shape, and mechanical properties such as elasticity cannot be evaluated.
The PWV method has been developed to measure theelasticity of the arterial wall.13–15) Propagation of the pressure
wave along the arterial tree is evaluated in this noninvasive
diagnostic method. Although it is useful for evaluation of
global elasticity such as from the heart to the femoral artery,regional elasticity cannot be evaluated.
In a series of studies of Kanai and coworkers,16–29) theelasticity of a local region has been evaluated by the ultrasonicphased-tracking method. In this method, a small change inthe thickness of the arterial wall due to the heartbeat was
accurately measured in each local region. From the measuredchanges in thickness and blood pressure, the regionalelasticity of the arterial wall was evaluated noninvasively.
In the measurement of the elasticity of the arterial wall by
this method, the boundaries of the arterial wall, between
which the elastic modulus is calculated, must be conﬁgured
beforehand. Currently, the media–adventitia and lumen–

intima boundaries of the anterior and posterior arterial walls
are manually determined by the operator. Since the manual
conﬁguration of the boundaries is time-consuming and results
in inter- and intra-observer variabilities, a fast automatedmethod is necessary for the detection of the boundaries.
In previous research, eﬀorts to detect the boundaries havebeen made by investigators worldwide. Touboul et al. in
1992 reported that IMT was automatically computed as the
distance between midpoints of the rising slopes of the twopulses from the lumen–intima and media–adventitia boundaries.30) In 1993, Gariepy et al. presented a system with which
the lumen–intima and media adventitia boundaries aredetected from changes in intensity based on default positions
that are traced by an operator.31) Selzer et al. in 1994
described a three-step method as follows: (1) approximate
echo boundaries are identiﬁed manually; (2) the approximateboundaries are used to guide the edge-ﬁnding algorithm;(3) the initial edges are tested in terms of edge strength, andweak edges are eliminated.32) Beux et al. also used thegradient of the intensity in the B-mode image of the artery.33)
The above methods all required a default setting by ahuman operator, which could bring variability owing tosubjective judgment. The above algorithms are based on
the evaluation of local features in a single image, either theintensity or the intensity gradient. These methods can detectboundaries precisely if a B-mode image is clear. Unfortunately, most B-mode images in a clinical situation containvarious types of noise such as speckle noise and multiple
reﬂection noise, making it diﬃcult to detect the boundaries
precisely. Several methods using a cost function have been
proposed for the detection of the boundaries, which arerobust against noise.34–40) Among these methods, Liang et al.
proposed a fully automated method using multiscale dynamic
programming (DP), in which the cost function is composedof echo intensity, echo intensity gradient, and change invertical distance.36)
In this paper, we propose an automated method for thedetection of arterial wall boundaries using multiscale DP for
automated measurement of elasticity, in which the correlation
between ultrasonic echoes in adjacent receive scan lines is
added to the cost function as a new feature. The correlation

07HF18-1© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.

10 mm

2 mm

1 sec

2 mm

[MPa]
4

2 mm

2
0

Zoom
ROI

(b)

(a)

Fig. 2. (Color online) Procedure for measuring the elasticity of the arterial
wall. (a) Zoomed B-mode image of carotid artery with manually traced
boundaries overlaid. (b) M-mode image in the red line in (a) with anterior
and posterior wall boundary tracking results overlaid. (c) Elasticity image.

Fig. 1. (Color online) Carotid artery B-mode image.

term enables more robust boundary detection to various typesof noise. Also, this method is used for the automateddetection of cardiac cycles, which is indispensable for the
automated measurement of elasticity. The eﬀectiveness of the
proposed method is validated using a phantom and applied tothe in vivo measurement of carotid arteries.
2.

Principles

2.1

Elasticity imaging method
To show the need for automated boundary detection, manual
trace processing for assessment of the elasticity of the carotid
arterial wall is shown in this section. The displacement ofan object is estimated from the quadrature demodulated
signal of the received ultrasonic RF echoes. Experiments
were conducted using ultrasonic diagnostic equipment,
FUJIFILM FAZONE M and an L10-5 linear type probe(with a frequency range from 5 to 10 MHz). The quadrature
demodulated signal was saved by a research data acquisition
tool named “IQscan”.41)
As shown in Fig. 1, a B-mode image of the carotid arteryis zoomed to obtain a high frame rate of 200 Hz for
acquisition of quadrature demodulated signals. The typical
M
C =2
X

(c)

zoom image size is 4 mm in width, 12 mm in height. The
sampling interval of the quadrature demodulated signal andthe pitch of scan lines were 0.035 and 0.12 mm, respectively.
For the zoomed image, the lumen–intima and media
adventitia boundaries of the anterior and posterior walls aremanually assigned. To reduce tracking error, an image frame
used for manual tracing is chosen close before the timing of
the minimum diameter. Figure 2(a) shows the zoomed Bmode image of the carotid artery with manually traced
boundaries overlaid.
With respect to Nl sampled points between the assigned
lumen–intima boundary and the media–adventitia boundary
along the l-th ultrasonic beam, the displacement of eachpoint, xi(t; l) (i = 1, 2, …, Nl ), is estimated using the phasedtracking method.16) For estimation of the displacement, using
the quadrature demodulated signal z(t; d; l) reﬂected at thetime t, at the depth d along the l-th ultrasonic beam, the phaseshift ΔΨi(t; l) between echoes in two consecutive frames isobtained from the complex cross-correlation function as
follows:

zðt þ Tr ; d þ xi ðt; lÞ þ mD; lÞ  z ðt; d þ xi ðt; lÞ þ mD; lÞm¼M =2

C;
exp½ j ^ i ðt; lÞ ¼  M =2

 X
C



zðt þ Tr ; d þ xi ðt; lÞ þ mD; lÞ  z ðt; d þ xi ðt; lÞ þ mD; lÞm¼M =2

ð1ÞC

where D and Tr are the interval of sampled points in the depth
direction and the pulse repetition interval, respectively, and +represents the complex conjugate. In Eq. (1), the number ofsampled points, Mc + 1, is set at 9 (= 0.41 µs) in consideration of the pulse length of 0.45 µs. In estimation of thephase shift by Eq. (1), the object position is tracked byintegration of the average velocity vi(t + Tr=2; l) during thepulse repetition interval Tr as follows:Trx^ i ðt þ Tr ; lÞ ¼ x^ i ðt; lÞ þ v^ i t þ ; l  Tr2
c0
¼ x^ i ðt; lÞ þ
 ^ i ðt; lÞ  Tr ;ð2Þ
2!0
where ω0 and c0 are the center angular frequency of theultrasonic pulse and the speed of sound, respectively.
Figure 2(b) shows the M-mode image in the red line on
Fig. 2(a) with tracking results for the anterior and posteriorwall boundaries overlaid.

From the displacement xi(t; l) estimated during therepetition interval Tr, which are set at the i-th point in thearterial wall along the l-th ultrasonic beam, the small change
in the thickness Δhi(t; l) of the arterial wall is obtained asfollows:
hi ðt; lÞ ¼ xi ðt; lÞ  xiþ1 ðt; lÞZt¼ ½vi ðt; lÞ  viþ1 ðt; lÞ dt:0

ð3Þ

From the ratio of the maximum decrease in thickness during
a cardiac cycle, Δhi,max(l) = maxt∣Δhi(t; l)∣, to the initial
thickness, h0(l), of the i-th layer, the maximum strain of thei-th layer is obtained by Δεi,max(l) = Δhi,max(l)=h0(l).
By assuming that the arterial wall is incompressible andthat the blood pressure is applied normal to each layer, theelastic modulus Eθ,i(l) of the i-th layer along the l-th beam isapproximately given by21)

07HF18-2© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.
Level 0(Original)

Gaussian pyramid to decompose to level 3 image

The coarsest scale (level 3 image) DP, detect lumen,anterior wall, and posterior wall

11.8 mm

Level 1

Coarse scale (expanded level 2 image) DP, detect media-adventitia boundary of anterior wall and posterior wall
Level 2Phased-tracking method, detect minimum diameter timing

Level 3

Fine scale (expanded level 1 image) DP, detect mediaadventitia boundary of anterior wall and posterior wall at
minimum diameter timing

3.8 mmFig. 4. Gaussian pyramid images.
Level 0 Expanded Expanded Expanded(Original) Level 1 Level 2 Level 3

Fine scale (expanded level 1 image) DP, detect lumen-intima boundary of anterior wall and posterior wall at
minimum diameter timing
Fig. 3. Outline of automated detection procedure.


3
2i0 ðlÞ
p1þ
;E;i ðlÞ ¼
8h0 ðlÞ "i;max ðlÞ

ð4Þ

where ρi0(l) and Δp are the initial inner radius of curvature of
the i-th layer along the l-th beam at a time t0 of the minimum
diameter and the pulse pressure (diﬀerence between thesystolic blood pressure and the diastolic blood pressure)measured at the upper arm, respectively. The image of theregional elasticity of the posterior wall is obtained, as shownin Fig. 2(c). We have mostly focused on the regionalelasticity of the posterior wall in our series of studies because
it is often diﬃcult to measure the regional elasticity of theanterior wall on clinical images because of multi-reﬂection
artifacts and sidelobe artifacts. As described above, thelumen–intima and media–adventitia boundaries of theanterior and posterior walls must be determined for measurement of the elasticity of the posterior wall by this method.
2.2

Automated artery boundary detection method
The multiscale DP algorithm36) is a technique for ﬁnding theboundaries by maximizing a certain cost function using
multiscale images. In the present study, by applying thistechnique to pre-scan-converted B-mode images, which areconstructed from the amplitude of the quadrature demodulated signal, the boundaries are detected. The reasons why
we use pre-scan-converted B-mode images are as follows:
introducing the correlation term between ultrasonic echoes inadjacent receive scan lines and reduction in computationalload owing to reduction in image size. In a pre-scanconverted B-mode image, the pixel spacing is 0.12 mm=pixelin the lateral direction and 0.035 mm=pixel in the axial
direction. An outline of the automated detection procedure isshown in Fig. 3.
2.2.1

11.8 mm

3.8 mmFig. 5. Levels of Gaussian pyramid expanded to size of original image.

deviation σ of 1 for normalized envelope signals is used forsmoothing.42) In the level 1 image, the standard deviation σof 1 and kernel size of 5 pixels in the axial direction werechosen to remove small noise and to detect the boundaries
with a high degree of accuracy and preserve the shape of an
ultrasound pulse, corresponding to approximately 10 pixels.
The kernel size in the lateral direction is longer than that in
the axial direction because a horizontally longer kernel makes
it easier to detect arterial wall boundaries, which are a
multilayer structure and often horizontally continuous. The
original image is repeatedly ﬁltered and subsampled togenerate the sequence of reduced resolution images. The
Gaussian pyramid images are shown in Fig. 4 and the levelsof the Gaussian pyramid expanded to the size of the original
image are shown in Fig. 5.
2.2.2 Cost function and dynamic programming In a gridof image M × N, consider a deformable polyline vector b
containing N nodes, one in each column composed of M
boundary candidates:b ¼ ðb1 ; b2 ; . . . ; bi1 ; bi ; . . . ; bN Þ;

Multiscale image Firstly, we use a Gaussian pyr-amid to decompose an ultrasound image into multiscale
images where a 5 × 5 Gaussian ﬁlter with the standard

ð5Þ

where bi−1 and bi are the horizontal neighboring points. At
point bi, a local cost C(bi) is deﬁned by

07HF18-3© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.
0.07

0.07

0.06

A

Erms [mm]Erms [mm]

0.06
0.05
0.04

B

A
0.05
0.04

B

0.03

0.03

C

C

(a)

(b)

0.02

0.02
0

0.125 0.25 0.375

0.5

0.625 0.75 0.875

01

0.125 0.25 0.375

0.625 0.75 0.875

1

w2

w1
0.07

0.07

0.06

0.06

A

Erms [mm]Erms [mm]

0.5

0.05
0.04

B

A
0.05
0.04

B

0.03

0.03

C

C

(c)

0.02

(d)

0.02
0

0.125 0.25 0.375

0.5

0.625 0.75 0.875

1

0

w3

0.125 0.25 0.375

0.5

0.625 0.75 0.875

1

w4

Fig. 6. (Color online) Averaged RMS error. (a) Projection of RMS error onto w1 axis. (b) Projection of RMS error onto w2 axis. (c) Projection of RMS error
onto w3 axis. (d) Projection of RMS error onto w4 axis.

Cðbi Þ ¼ w1 Iðbi Þ þ w2 Gðbi Þ þ w3 Rðbi1 ; bi Þþ w4 Dðjbi1  bi jÞ;

ð6Þ

where I, G, R, and D are the cost function terms, which
represent respective boundary features, and w1, w2, w3, andw4 are their weight factors, which reﬂect the relative
importance of each boundary feature in each image level.
Term I represents the echo intensity feature, which is anormalized diﬀerential value between the sum intensity of
pixels below the interface and the sum intensity of pixels
above the interface. The term G represents the echo intensity
gradient feature, which is the normalized value of thedownward intensity slope. The term R represents the lateralsimilarity to the feature in shape around the boundary. In
the present study, R is employed as the cross correlationcoeﬃcient of ultrasonic signals around bi−1 and bi. The crosscorrelation window length is almost equal to the pulse length.
The term D is the vertical position continuity of theboundary, which prefers a smoother line,
8< jbi1  bi j
for jbi1  bi j ≦ LDðjbi1  bi jÞ ¼
; ð7Þ
L
:
1
for jbi1  bi j > L
where L is the supposed maximum change in the verticaldistance between the adjacent lines and is set to 5 pixels(0.175 mm) here to detect a relatively sharp change in thevertical distance of lesions such as plaques. The relationshipamong the terms I, G, R, and D and the eﬀect of the term R
are explained in Sect. 2.2.4.
The cost function is deﬁned as the sum of local costs along
the polyline.
N
X
Csum ¼Cðbi Þ:ð8Þ
i¼1

The optimum boundary is searched to maximize the costfunction Csum. An exhaustive whole space M × N searchneeds a computational complexity O(MN), which requires
numerous calculations and could need a longer time thanmanual tracing. To reduce calculation time signiﬁcantly,
we search the optimum boundary by applying the DP
procedure, and the computational complexity is reduced toO(N × M2).35,43,44)
Determining the weight factors In this study, 20
subject images were analyzed. One-half of the images wereused as training data for determining the weight factors inthis section, the other half of the images were analyzed in
Sect. 3.2. The cost function for detecting one interface
includes four terms. Because we use the constraint w1 +w2 + w3 + w4 = 1, we only determine three of them. The
combinations of ðw1 ; w2 ; w3 Þ, in steps of 1=8, such thatw1 + w2 + w3 ≤ 1, have 165 patterns. For the values ofðw1 ; w2 ; w3 Þ, we apply the DP procedure to the training
data and compute the root mean square (RMS) error between
the automatically and manually assigned boundaries by aprofessional expert. The RMS error eRMS is deﬁned by thedepth dauto(l) and dman(l) of the automatically detectedboundary and manually assigned boundary at the beam
position l as follows:
vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
u
Nb
u1 X
ð9Þ
eRMS ¼ t
jdauto ðlÞ  dman ðlÞj2 ;Nb l¼1
2.2.3

where Nb is the number of beam positions. The values ofðw1 ; w2 ; w3 ; w4 Þ associated with the lowest RMS error foreach boundary are determined as the weight factors. As an
example, we present the relationships between average RMS
error among the 10 training images and the values ofðw1 ; w2 ; w3 ; w4 Þ for the media–adventitia boundary of theposterior wall of the expanded level 1 image in Fig. 6.

07HF18-4© 2015 The Japan Society of Applied Physics

Y. Miyachi et al.

0.07

0.07

0.06

0.06
Erms [mm]Erms [mm]Jpn. J. Appl. Phys. 54, 07HF18 (2015)

0.05
0.04

0.05
0.04
0.03

0.03

(b)

(a)
0.02

0.02
0

0.125 0.25 0.375

0.5

0.625 0.75 0.875

01

0.125 0.25 0.375

0.625 0.75 0.875

1

w2

w1
0.07

0.07

0.06

0.06
Erms [mm]Erms [mm]

0.5

0.05
0.04

0.05
0.04
0.03

0.03

(d)

(c)
0.02

0.02
0

0.125 0.25 0.375

0.5

0.625 0.75 0.875

01

0.125 0.25 0.375

0.5

0.625 0.75 0.875

1

w4

w3

0.07

0.07

0.06

0.06

Erms [mm]Erms [mm]Fig. 7. (Color online) RMS error for image without speckle noise, in which boundaries are clearly seen. Data with minimum RMS errors are indicated by
red diamonds. (a) Projection of RMS error onto w1 axis. (b) Projection of RMS error onto w2 axis. (c) Projection of RMS error onto w3 axis. (d) Projection of
RMS error onto w4 axis.

0.05
0.04

0.05
0.04
0.03

0.03

(b)

(a)

0.02

0.02
0

0.125 0.25 0.375

0.5

0.625 0.75 0.875

01

0.125 0.25 0.375

0.625 0.75 0.875

1

w2

w1
0.07

0.07

0.06

0.06

Erms [mm]Erms [mm]

0.5

0.05
0.04
0.03

0.05
0.04
0.03

(c)

(d)

0.02

0.02
0

0.125 0.25 0.375

0.5

0.625 0.75 0.875

1

0

w3

0.125 0.25 0.375

0.5

0.625 0.75 0.875

1

w4

Fig. 8. (Color online) RMS error for image with speckle noise, in which boundaries are partly unclear. Data with minimum RMS errors are indicated by red
diamonds. (a) Projection of RMS error onto w1 axis. (b) Projection of RMS error onto w2 axis. (c) Projection of RMS error onto w3 axis. (d) Projection of RMS
error onto w4 axis.

Figures 6(a)–6(d) show the projection of eRMS onto the w1,
w2, w3, and w4 axes respectively. The eRMS axis is zoomedaround the minimum values. eRMS is minimum when w1, w2,
w3, and w4 are non-zero values. This suggests that it is betterto consider all the terms I, G, D, and R. The values ofðw1 ; w2 ; w3 ; w4 Þ associated with the lowest RMS error were(0.375, 0.125, 0.375, 0.125) and these weight factors wereused for the media–adventitia boundary of the posterior wallin the analysis in Sect. 3.2.
To investigate the tendencies of eRMS and the weightfactors, the projection results of eRMS onto w1, w2, w3, andw4 of two subjects are shown in Figs. 7 and 8. Figure 7 is fora clear boundary image without speckle noise and Fig. 8 isfor a partly unclear boundary image with speckle noise.

07HF18-5© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)
3.8 mm(33 lines)

w1 I (bi )Intensity
I (bi )

4.4 mm(126
sampledpoints)Y. Miyachi et al.

+ w2 G (bi )Valley ofSpeckle Noise

0
y

x

0

y
(sampling)y(sampling)x (Line)x (Line)

(a)

(b)

(a)

w1 I (bi )+ w2 G (bi )Fig. 9. (Color online) (a) Zoomed level 1 image of posterior wall. (b) 3D
plot of level 1 image intensity around posterior media–adventitia boundary.

3.8 mmy(sampling)

4.4 mmx (Line)(b)
w1 I (bi )

(a)

(b)

(c)

+ w2 G (bi )(d)Fig. 10. (Color online) Manually traced and automatically detectedboundary. (a) Boundary manually traced by a professional expert.
(b) Automatically detected boundary considering only terms I and G
calculated using values of point A. (c) Automatically detected boundary
considering terms I, G, and D calculated using values of point B.
(d) Boundary automatically detected by our method considering terms I, G,
D, and R calculated by using values of point C.

y(sampling)

In both ﬁgures, eRMS is minimum at weight factors of(0.375, 0.125, 0.375, 0.125) plotted with red diamonds. The
dependence of weight factors observed in Fig. 8 is stronger
than that observed in Fig. 7.
Effect of the terms I, G, D, and R So far, we havedescribed the DP method. Here, we explain the eﬀect of theterms I, G, D, and R using Fig. 6. Point A is (0.75, 0.25, 0, 0)of ðw1 ; w2 ; w3 ; w4 Þ, which leads to the minimum eRMS when
w3 and w4 are zero. Point B is (0.5, 0.25, 0, 0.25) ofðw1 ; w2 ; w3 ; w4 Þ, which leads to the minimum eRMS when
w3 is zero. Point C is (0.375, 0.125, 0.375, 0.125) ofðw1 ; w2 ; w3 ; w4 Þ, which leads to the minimum eRMS in all
combinations of ðw1 ; w2 ; w3 ; w4 Þ.
We explain the reason why the values of point C give theminimum eRMS using a typical image, where the detection
result using values of point C leads to a much lower eRMS
than that using values of points A and B. Figures 9(a) and9(b) show the zoomed level 1 image of the posterior wall andthe three-dimensional (3D) plot of level 1 image intensity
around the posterior media–adventitia boundary. As shown in
Figs. 9(a) and 9(b), the valley of speckle noise exists aroundthe left half 1–3 lines and 10–16 lines.
Figure 10(a) shows the boundary manually traced by aprofessional expert. We can understand that the professionalx (Line)(c)Fig. 11. (Color online) 3D plot of w1I + w2G. (a) 3D plot of w1I + w2G
calculated using values of point A. (b) 3D plot of w1I + w2G calculatedusing values of point B. (c) 3D plot of w1I + w2G calculated using values ofpoint C.

2.2.4

expert traces the boundary such that speckle noise edges
are avoided and boundary continuity is given importance.
Figure 10(b) shows the boundary automatically detectedusing the values of point A and Fig. 11(a) shows the 3D plotof w1I + w2G. The use of the values of point A detects
speckle noise edges by mistake and leads to a higher eRMS,
0.105 mm, because the intensity and gradient at the boundarywith speckle noise are higher than those at the manually
traced boundary. Because the values of point A, which are
w3 and w4, are zero, the continuity of the boundary is notconsidered. Figure 10(c) shows the automatically detectedboundary using values of point B, and Fig. 11(b) shows the3D plot of w1I + w2G. The use of the values of point B also
detects speckle noise edges by mistake and leads to a higher
eRMS, 0.135 mm. This is because the manually traced
boundary indicates lower values of intensity and gradientthan the speckle noise edges at the lines with speckle noise,and the speckle noise edge depth changes continuously,
which means the change in the speckle noise edge depth

07HF18-6© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.

3.8 mm

w1 I (bi )+ w2 G (bi )4.4 mm

(a)

(b)

(c)

y
(sampling)(d)Fig. 12. (Color online) Manually traced and automatically detectedboundaries for another subject. (a) Boundary manually traced by aprofessional expert. (b) Automatically detected boundary considering only
terms I and G calculated using values of point A. (c) Automatically detectedboundary considering terms I, G, and D calculated using values of point B.
(d) Boundary automatically detected by our method considering terms I, G,
D, and R calculated using values of point C.

between adjacent lines is small. Figure 10(d) shows theboundary automatically detected using the values of point C,
and Fig. 11(c) shows the 3D plot of w1I + w2G. The use of
the values of point C leads to a lower eRMS, 0.028 mm, which
is close to the manually traced result. As one more example,
the manually traced and detection results of another subject
are shown in Figs. 12(a)–12(d) and the 3D plots of w1I +w2G are shown in Figs. 13(a)–13(c).
These results show that the term R, which represents thelateral similarity between the waveforms of ultrasonic echoes,
is eﬀective even in the boundary with speckle noise. Thisis because the lateral similarity between adjacent lines isconserved even in the boundary with speckle noise, as we cansee in 1–3 lines and 10–16 lines in Fig. 9(b), although theintensity and gradient are relatively lower in the boundarywith speckle noise. The results of detection using the valuesof points B and C indicate that a professional expert traces
the boundary in terms of not only the continuity of thechange in the vertical distance but also the continuity of the
waveform especially in the boundary with speckle noise.
Regardless of whether the term R is used, detection results
are similar at clear boundary lines without speckle noise suchas 16–33 lines, because the intensity and gradient are higher.
However, at boundary lines with speckle noise such as 10–16
lines, the intensity and gradient are lower and the speckle
noise edge depth changes continuously; only the detection
considering the term R shows results close to those in which
the boundary is manually assigned.
Detection of the minimum diameter timing To detect the minimum diameter timing at which precise boundary
detection is conducted, ﬁrstly, we detect the lumen, anteriorwall, and posterior wall roughly using the coarsest scale(level 3) image. In the level 3 image, the center depth of thelumen is easily determined as the depth where the mean
intensity obtained by averaging the intensity proﬁle laterally
is minimum. The level 3 image is divided into the upper partand lower part at the depth of the centered lumen, and DP is
conducted in the upper part to detect the anterior wall and in
the lower part to detect the posterior wall. Figure 14(a) shows
the level 3 image, and Fig. 14(b) shows the level 3 image
with the automatically detected lumen, anterior, and posteriorx (Line)(a)w1 I (bi )+ w2 G (bi )y(sampling)x (Line)(b)
w1 I (bi )+ w2 G (bi )y(sampling)x (Line)(c)Fig. 13. (Color online) 3D plot of w1I + w2G for another subject. (a) 3D
plot of w1I + w2G calculated using values of point A. (b) 3D plot of
w1I + w2G calculated using values of point B. (c) 3D plot of w1I + w2G
calculated using values of point C.

3.8 mmAnterior
Wall

2.2.5

Lumen

11.8 mm

Posterior
Wall

(a)

(b)

(c)

(d)

Fig. 14. (Color online) (a) Level 3 image. (b) Level 3 image with
automatically detected lumen, anterior, and posterior wall boundaries
overlaid. (c) Extended level 2 image. (d) Extended level 2 image with
automatically detected anterior media–adventitia boundary and posterior–
media adventitia boundary overlaid.

07HF18-7© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.
3.8 mm2 mmAnterior media-adventitia boundary
Anterior lumen-intima boundary

11.8 mmchange 0.3
in diameter[mm]
0

Posterior lumen-intima boundary
Posterior media-adventitia boundary

(b)

(a)

(a)

(c)
1

2
time [s]3Fig. 15. (Color online) (a) Extended level 2 image with automaticallydetected anterior media–adventitia boundary and posterior media–adventitia
boundary overlaid. (b) M-mode image in the red line in Fig. 11(a) withanterior and posterior wall media–adventitia boundary tracking results
overlaid. (c) Change in vascular diameter.

wall boundaries overlaid. Secondly, the media–adventitia
boundary is detected because echoes from the media–
adventitia boundary are larger and easier to detect than thosefrom the lumen–intima boundary. The anterior media–
adventitia boundary is detected in the expanded level 2
image by DP searching around the depth of the arterial wall
detected in the level 3 image, and the posterior media–
adventitia boundary is detected in the expanded level 2
image by DP searching around the depth of the posterior wall
detected in the level 3 image. The DP searching range iswithin ±1 mm from the depth of the detected boundary inlevel 3, because a typical gap between the coarsely and ﬁnely
determined boundary positions is less than 0.35 mm of theultrasonic pulse length. An expanded level 2 image is shownin Fig. 14(c) and an expanded level 2 image with theautomatically detected anterior media–adventitia boundaryand posterior media–adventitia boundary overlaid is shown in
Fig. 14(d). Lastly, to detect the change in diameter and local
minimum diameter timing, the phased-tracking method isapplied to the quadrature demodulated signal from the depthsof the detected anterior media–adventitia boundary andposterior media–adventitia boundary in the expanded level2 image. Figure 15(a) shows the level 2 image with theautomatically detected anterior media–adventitia boundaryand posterior media adventitia boundary overlaid, andFig. 15(b) shows the M-mode image in the red line in
Fig. 15(a) with the anterior and posterior wall media–
adventitia boundaries tracking results overlaid. The change
in diameter estimated by the phased-tracking method isshown in Fig. 15(c). The local minimum diameter timings areindicated by the green lines in Figs. 15(b) and 15(c).
2.2.6

Boundary detection at minimum diameter timing
Using the expanded level 1 image at local minimum diameter
timing, ﬁne-scale DP detection of the anterior media–
adventitia boundary, anterior lumen–intima boundary, posterior media–adventitia boundary, and posterior lumen–intima
boundary is conducted. Firstly, the media–adventitia boundary is detected on the basis of the boundary positions
detected in a coarser scale (expanded level 2) image.

(b)

Fig. 16. (Color online) (a) Expanded level 1 image at minimum diameter
timing. (b) Expanded level 1 image with automatically detected anterior
lumen–intima boundary, anterior media–adventitia boundary, posterior
lumen–intima boundary, and posterior media–adventitia boundary at
minimum diameter timing.

transducer

beef
silicon rubber
CCD
tube
Laser
micrometerwater tankstepping
motor
pump

Fig. 17. (Color online) Arterial-wall-mimicking-phantom made ofsilicone rubber tube pulsated by stepping motor pump.

Subsequently, the lumen–intima boundary is detected bysearching 2.5 mm inside the media–adventitia boundary
position detected. The search range of 2.5 mm is chosen to
cover a typical intima–media thickness, which is less than
2 mm. The expanded level 1 image is shown in Fig. 16(a),and the expanded level 1 image with the automaticallydetected anterior lumen–intima boundary, anterior media–
adventitia boundary, posterior lumen–intima boundary, andposterior media–adventitia boundary is shown in Fig. 16(b).
3.

Results

3.1

Basic experiments using silicone tube phantom
To evaluate the accuracy of the detected minimum diameter
timing and boundaries, we prepared a phantom mimicking
the arterial wall made of a silicone rubber tube, which was
pulsated by a stepping motor pump, as shown in Fig. 17. The
thickness of the silicone rubber tube is 1.1 mm and the innerdiameter is 8 mm, which are similar to those of the carotid
artery. The material of the silicone rubber tube includes
graphite powder (weight ratio, 5%), acting as ultrasoundscatters. To evaluate the change in external diameter by amethod other than the proposed method, a CCD laser
micrometer (Keyence IG-028) is located across a transparent
water tank. To mimic random speckle noise in the B-mode
used in clinical settings, a thick piece of beef was located on
the silicone rubber tube.

07HF18-8© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.
2 mm

2 mm

(a)

2 mm

2 mm

Min. diameter timing
detected by phasedtracking methodchange 10.4
in diameter[mm]
10change 10.4
in diameter
[mm]
10
0

(b)

(d)

(a)
(c)
0.6
time [s]1.2

Fig. 18. (Color online) (a) Arterial-wall-mimicking-phantom M-mode
image with anterior and posterior wall media–adventitia boundary tracking
results overlaid. (b) Change in diameter measured by phased-tracking
method. (c) Change in diameter measured using CCD laser micrometer.
(d) Arterial-wall-mimicking-phantom B-mode image with automaticallydetected anterior lumen–intima boundary, anterior media–adventitia
boundary, posterior lumen–intima boundary, and posterior media adventitia
boundary.

Figures 18(a)–18(c) show M-mode images with the
tracking results with respect to the anterior and posteriorwall media–adventitia boundaries overlaid, the change indiameter measured by the phased-tracking method, and thechange in diameter measured by the CCD laser micrometer,
respectively. The diameter was calculated by the phasedtracking method considering a sound speed of 1030 m=s forthe silicone rubber tube and that of 1480 m=s for water.
Figure 18(d) shows a scan-converted B-mode image withthe automatically detected anterior lumen–intima boundary,anterior media–adventitia boundary, posterior lumen–intima
boundary, and posterior media adventitia boundary. Comparing Fig. 18(b) with Fig. 18(c), the diﬀerence in minimum
diameter timings between the phased-tracking method andthe method using a CCD laser micrometer is less than
0.015 s. Figures 19(a)–19(d) show a scan-converted B-mode
image, a scan-converted B-mode image with automaticallydetected boundaries considering only terms I and G overlaid,a scan-converted B-mode image with the automaticallydetected boundaries considering terms I, G, and D overlaid,
and a scan-converted B-mode image with the automaticallydetected boundaries considering terms I, G, D, and R overlaid
at the minimum diameter timing, respectively. The averageand maximum thicknesses of the silicone rubber tube arecalculated considering a sound of 1030 m=s for the silicone
rubber tube. The average thicknesses of the anterior wall are
1.18 mm in Fig. 19(b), 1.09 mm in Fig. 19(c), and 1.09 mm
in Fig. 19(d). Also, the average thicknesses of the posteriorwall are 1.05 mm in Fig. 19(b), 1.04 mm in Fig. 19(c), and1.06 mm in Fig. 19(d). The maximum thicknesses of theanterior wall are 1.44 mm in Fig. 19(b), 1.37 mm in

(b)

(c)

(d)

Fig. 19. (Color online) (a) Scan-converted B-mode image of arterial-wallmimicking phantom at minimum diameter timing. (b) Scan-converted Bmode image with automatically detected boundaries considering only terms
I and G overlaid. (c) Scan-converted B-mode image with automaticallydetected boundaries considering terms I, G, and D overlaid. (d) Scanconverted B-mode image with automatically detected boundaries consideringterms I, G, D, and R overlaid.

Fig. 19(c), and 1.12 mm in Fig. 19(d). The maximum
thicknesses of the posterior wall are 1.14 mm in Fig. 19(b),1.14 mm in Fig. 19(c), and 1.14 mm in Fig. 19(d). The
thickness of the silicone rubber tube is measured to be
1.1 mm under a microscope and the surface is conﬁrmed to
be ﬂat; the estimated maximum thicknesses of the anteriorwall in Figs. 19(b) and 19(c) have a margin of error. On theother hand, the estimated average and maximum thicknesses
in Fig. 19(d) are well in agreement with the results ofmeasurement under a microscope.
3.2

In vivo experimental results for human carotid artery
The proposed method was applied to the in vivo measurement of carotid arteries, and the root mean square errors
between the proposed method, and the manual assignmentof the boundaries were evaluated in 10 subjects with andwithout atherosclerotic plaques. Examples of images of a
subject without atherosclerotic plaques are shown in Fig. 20
and those of a subject with atherosclerotic plaques are shownin Fig. 21. In both ﬁgures, (a), (b), (c), (d), and (e) show ascan-converted B-mode image, a scan-converted B-mode
image with the automatically detected boundaries considering only terms I and G overlaid, a scan-converted B-mode
image with the automatically detected boundaries considering terms I, G, and D overlaid, a scan-converted B-mode
image with the automatically detected boundaries considering terms I, G, D, and R overlaid, and a scan-converted
B-mode image with the boundaries manually traced by aprofessional expert at the minimum diameter timing,
respectively. In Fig. 20, the relatively clear boundaries,
anterior media–adventitia boundary, anterior lumen–intima
boundary, and posterior lumen–intima boundary, are very
close to the manually assigned boundaries for all automateddetection methods. On the other hand, for the posterior
media–adventitia boundary, in which speckle noise ispresent, only the proposed method is in good agreement

07HF18-9© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.
2 mm
2 mm

(a)

(b)

(c)

(d)

(e)

Fig. 20. (Color online) Examples of images of subject without plaques. (a) Scan-converted B-mode image. (b) Scan-converted B-mode image with
automatically detected boundaries considering only terms I and G overlaid. (c) Scan-converted B-mode image with automatically detected boundaries
considering terms I, G, and D overlaid. (d) Scan-converted B-mode image with automatically detected boundaries considering terms I, G, D, and R overlaid.
(e) Scan-converted B-mode image with the boundaries manually traced by a professional expert.

2 mm
2 mm

(a)

(b)

(c)

(d)

(e)

Fig. 21. (Color online) Examples of images of a subject with plaques. (a) Scan-converted B-mode image. (b) Scan-converted B-mode image with
automatically detected boundaries considering only terms I and G overlaid. (c) Scan-converted B-mode image with automatically detected boundaries
considering terms I, G, and D overlaid. (d) Scan-converted B-mode image with automatically detected boundaries considering terms I, G, D, and R overlaid.
(e) Scan-converted B-mode image with the boundaries manually traced by a professional expert.

with the manually traced boundaries. In Fig. 21, the proposed
method shows a smoother plaque boundary, which is closer
to the manual tracing than the other methods. In 10 subjects,
the root mean square errors between all methods and themanual assignment are 0.069 mm for the method considering
only terms I and G, 0.050 mm for the method consideringterms I, G, and D, and 0.034 mm for the proposed method
considering terms I, G, D, and R. Compared with the othermethods, the proposed method indicates smaller root meansquare errors.
4.

Discussion

In conventional arterial wall boundary detection procedures,
intensity and its gradient terms of each scan line are used.

Although such simple methods lead to precise detection if a
B-mode image is clear, they are not suﬃcient for most Bmode images in clinical settings because they contain various
types of noise such as speckle noise. Therefore, several
methods using cost function have been proposed for robust
noise detection of the boundaries. Liang et al. proposed afully automated method using multiscale DP, in which thecost function is composed of echo intensity, echo intensity
gradient, and change in vertical distance.36) However, thismethod is also insuﬃcient to balance strong speckle noises
with abrupt boundary changes in the vertical distance in casessuch as plaques. As shown in Sect. 3, adding the viewpointof lateral similarity feature with respect to waveforms ofultrasonic echoes to cost function enables the balance

07HF18-10© 2015 The Japan Society of Applied Physics

Jpn. J. Appl. Phys. 54, 07HF18 (2015)

Y. Miyachi et al.

between noise robustness and detection of the abrupt change
in the boundary. It is important and eﬀective to consider not
only the continuity in the vertical position but also the lateralcontinuity with respect to waveforms of ultrasonic echoes forthe boundary detection.
5.

Conclusions

We have proposed an automated arterial wall boundary
detection method using multiscale DP, in which the costfunction includes the correlation term between ultrasonicechoes in adjacent receive scan lines. The accuracy of the
proposed method was validated using a phantom mimicking
the arterial wall made of a silicone rubber tube pulsated by astepping motor pump. The proposed method was also applied
to the in vivo measurement of carotid arteries, and the results
were compared with results manually traced by a humanexpert. The proposed method was much closer to the results
obtained by a human expert than to the results of the simple
method using only intensity and its gradient in each scan line
and the DP method using only the change in vertical distancein successive scan line. The correlation term enables
boundary detection to be more robust against noises. The
proposed method shows potential for detecting arterial wallboundaries automatically and may be useful for elasticity
imaging.

1) H. Tanaka, M. Nishino, M. Ishida, R. Fukunaga, and K. J. Sueyoshi, Stroke
23, 946 (1992).
2) R. T. Lee, A. J. Grodzinsky, E. H. Frank, R. D. Kamm, and F. J. Schoen,
Circulation 83, 1764 (1991).
3) H. M. Loree, A. J. Grodzinsky, S. Y. Park, L. J. Gibson, and R. T. Lee,
J. Biomech. 27, 195 (1994).
4) P. C. G. Simons, A. Algra, M. L. Bots, D. E. Grobbee, and Y.
van der Graaf, Circulation 100, 951 (1999).
5) J. T. Salonen and R. Salonen, Circulation 87, 1156 (1993).
6) P. Pignoli, E. Tremoli, A. Poli, P. Oreste, and R. Paoletti, Circulation 74,
1399 (1986).
7) E. de Groot, G. K. Hovingh, A. Wiegman, P. Duriez, A. J. Smit, J.-C.
Fruchart, and J. J. P. Kastelein, Circulation 109, III-33 (2004).
8) S. Kageyama, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 52, 07HF04(2013).
9) C. Arihara, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 45, 4727
(2006).
10) T. Kaneko, H. Hasegawa, and H. Kanaii, Jpn. J. Appl. Phys. 46, 4881
(2007).
11) K. Kitamura, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 51, 07GF08(2012).
12) Y. Nagai, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 53, 07KF19(2014).
13) J. Blacher, R. Asmar, S. Djane, G. M. London, and M. E. Safar,

Hypertension 33, 1111 (1999).
14) P. Hallock, Arch. Int. Med. 54, 770 (1934).
15) M. Saito, Y. Yamamoto, Y. Shibayama, M. Matsukawa, Y. Watanabe, M.
Furuya, and T. Asada, Jpn. J. Appl. Phys. 50, 07HF10 (2011).
16) H. Kanai, H. Hasegawa, N. Chubachi, Y. Koiwa, and M. Tanaka, IEEE
Trans. Ultrason. Ferroelectr. Freq. Control 43, 791 (1996).
17) H. Hasegawa, H. Kanai, N. Hoshimiya, N. Chubachi, and Y. Koiwa, Jpn. J.
Appl. Phys. 37, 3101 (1998).
18) H. Hasegawa, H. Kanai, N. Hoshimiya, and Y. Koiwa, Jpn. J. Appl. Phys.
39, 3257 (2000).
19) H. Hasegawa, H. Kanai, and Y. Koiwa, Jpn. J. Appl. Phys. 41, 3563 (2002).
20) H. Hasegawa, H. Kanai, N. Chubachi, and Y. Koiwa, J. Acoust. Soc. Jpn.
53, 346 (1997) [in Japanese].
21) M. Cinthio, H. Hasegawa, and H. Kanai, Proc. IEEE Int. Ultrasonics Symp.,
2007, p. 997.
22) H. Hasegawa and H. Kanai, Jpn. J. Appl. Phys. 43, 3197 (2004).
23) J. Inagaki, H. Hasegawa, H. Kanai, M. Ichiki, and F. Tezuka, Jpn. J. Appl.
Phys. 44, 4593 (2005).
24) J. Tang, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 44, 4588 (2005).
25) K. Ikeshita, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 47, 4165
(2008).
26) K. Ikeshita, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 48, 07GJ10(2009).
27) K. Ikeshita, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 50, 07HF08(2011).
28) K. Ikeshita, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 51, 07GF14(2012).
29) M. Sato, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 53, 07KF03(2014).
30) P.-J. Touboul, P. Prati, P.-Y. Scarabin, V. Adrai, E. Thibout, and P.
Ducimetière, J. Hypertension 10, 37 (1992).
31) J. Gariepy, M. Massonneau, J. Levenson, D. Heudes, and A. Simon,
Hypertension 22, 111 (1993).
32) R. H. Selzer, H. N. Hodis, H. Kwong-Fu, W. J. Mack, P. L. Lee, C. R. Liu,and C. H. Liu, Atherosclerosis 111, 1 (1994).
33) F. Beux, S. Carmassi, M. V. Salvetti, L. Ghiadoni, Y. Huang, S. Taddei,
and A. Salvetti, Ultrasound Med. Biol. 27, 1621 (2001).
34) H. Hasegawa, H. Kanai, and Y. Koiwa, IEEE Trans. Ultrason. Ferroelectr.
Freq. Control 51, 93 (2004).
35) T. Gustavsson, Q. Liang, I. Wendelhag, and J. Wikstrand, Proc. IEEE
Computers Cardiology, 1994, p. 297.
36) Q. Liang, I. Wendelhag, J. Wikstrand, and T. Gustavsson, IEEE Trans.
Med. Imaging 19, 127 (2000).
37) I. Wendelhag, Q. Liang, T. Gustavsson, and J. Wikstrand, Stroke 28, 2195
(1997).
38) R. Rocha, A. Campilho, J. Silva, E. Azevedo, and R. Santos, Image Vision
Comput. 28, 614 (2010).
39) N. Ibrahim, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 51, 07GF07(2012).
40) N. Ibrahim, H. Hasegawa, and H. Kanai, Jpn. J. Appl. Phys. 52, 07HF03(2013).
41) L. Y. L. Mo, D. DeBusschere, W. Bai, D. Napolitano, A. Irish, S.
Marschall, G. W. McLaughlin, Z. Yang, P. L. Carson, and J. B. Fowlkes,
Proc. IEEE Int. Ultrasonics Symp., 2007, p. 2259.
42) P. J. Burt and E. H. Adeloson, IEEE Trans. Commun. 31, 532 (1983).
43) R. Bellman, Dynamic Programming (Princeton University Press, Princeton,
NJ, 1962) p. 180.
44) H. Satoh, H. Kanai, and N. Chubachi, Nihon Onkyo Gakkaishi 50, 11(1994) [in Japanese].

07HF18-11© 2015 The Japan Society of Applied Physics