Kluwer - Handbook of Biomedical Image Analysis Vol

232	Frangi, Laclaustra, and Yang

Figure 5.1: Experimental setup.

Each sequence has about 1200 frames and a duration of around 20 min, acquiring each second the last telediastolic frame previously hold. This provides a fixed sampling rate irrespective of heart rate, which means a substantial benefit for clinical interpretation, as different stimulus are applied on a time basis along the clinical test. As the dynamics of endothelial vasodilation is much slower than changes happening between cardiac cycles, with this sampling rate, missing information from one heartbeat or using one heartbeat twice does not affect, in practice, the results.

FMD is the vasodilation response to hyperemia after a transitory distal ischemia induced in the forearm using a pneumatic cuff distal to the probe (Fig. 5.1). The dilation mediated by a chemical vasodilator, the nitroglycerin, or nitroglycerin-mediated dilation (NMD), is also registered. Accordingly, five phases of the medical test can be distinguished in each sequence (see Fig. 5.2).

Rest baseline (B1). Initial rest state preceding distal ischemia. Presents the best image quality in the whole sequence and lasts for about 1 min.

Distal ischemia (DI). The cuff is inflated and, therefore, the image quality is usually the worst in the sequence. It takes approximately 5 min. This phase ends when the cuff pressure is released.

Computerized Analysis and Vasodilation Parameterization of FMD

233

Figure 5.2: A whole typical examination can be divided into several segments along the sequence: rest baseline (B1), distal ischemia (DI), flow-mediated dilation (FMD), post-FMD baseline (B2), and nitroglycerin-mediated dilation (NMD).

Flow-mediated dilation (FMD). Response to reactive hyperemia. The maximum %FMD of healthy subjects has a mean value around 5%, and it takes place about 60 sec after the cuff is released [5]. At the beginning of this phase, and coincident with the release of the cuff, it is common to have important motion artifacts. The duration of this phase depends on each patient.

Post-FMD baseline (B2). The artery diameter returns to a steady state. This state does not necessary have to coincide with that of B1.

Nitroglycerin-mediated dilation (NMD). Response to the sub lingual supply of a fixed dose of this chemical vasodilator, which is made 8 min later than the cuff release.

5.3Registration-Based FMD Estimation

5.3.1 Algorithm Overview

Our technique assumes that the vasodilation that takes place between two frames can be modeled by a constant scaling in the direction normal to the artery. This scale factor is obtained by means of image registration.

234	Frangi, Laclaustra, and Yang

Figure 5.3: Overview of the proposed two-stage method: after motion compensation is carried out by recovering a rigid motion model, the vasodilation is measured by computing the scaling factor along the normal to the artery that best matches the two analyzed frames.

A reference frame is selected from the beginning of the sequence. All the other frames are registered to this reference frame. Changes in the relative position between the patient and the transducer are quite common during a whole examination, which may take up to 20 min. To elude wrong anatomical correspondences, motion compensation becomes necessary, and a rigid image registration technique is used to this end.

Structures surrounding the artery in the image may be important to resolve potential ambiguities in the longitudinal alignment between two frames, which occur because of the lineal nature of the arterial walls. On the other hand, extraluminal structures introduce artifacts when measuring the vasodilation since they do not necessarily deform in the same way as the artery does. Therefore, they should be taken into account when retrieving the global rigid motion information of the model, while arterial vasodilation estimation should only consider the artery deformation.

Our technique proceeds in two phases as summarized in Fig. 5.3: motion compensation and dilation assessment. The first phase uses the original frames and rigid image registration to recover a rigid motion model. Translation and rotation parameters are used to initialize the subsequent phase of vasodilation estimation. This second stage employs an affine registration model. To avoid artifacts when measuring arterial vasodilation it is convenient to remove background extra luminal structures by padding them out from the reference frame. Preprocessing of this frame also requires repositioning it so that the artery is normal to the scaling direction, that is to say, aligned with the horizontal axis, since our model searches for a vertical scaling factor (see Fig. 5.4). Both operations are performed manually on the reference frame. Manual masking only requires to roughly draw two lines in the reference frame and repositioning, to align a line with the direction of the artery, a process that is simple and takes only a few seconds per image sequence.

Computerized Analysis and Vasodilation Parameterization of FMD

235

Figure 5.4: Preprocessing applied to the reference frame before the phase of vasodilation assessment. (a) Original reference frame and (b) the reference frame after alignment to the horizontal axis and padding out of background structures.

Temporal continuity is enforced in both phases by means of recursive filters in the registration parameter space prior to registering each new frame.

In the next two subsections the registration algorithm and its initialization, enforcing temporal continuity, are discussed.

5.3.2 Registration Algorithm

5.3.2.1 Motion and Vasodilation Models

Registering image B onto image A requires finding a transformation T (·) that maps B into A by maximizing a registration measure M(·) as indicated in Eq. (5.1). The similarity measure is computed over all points, P, of the overlap region of both images

M( A( P), B(T ( P)))

(5.1)

Our motion model between the original (x, y) and transformed (x , y ) coordinates is a rigid transformation of the form

y		=	sin θ	cos θ	·	y	+	ty
x			cos θ	−sinθ		x		tx	(5.2)

As we are imaging only a small and roughly straight vessel segment, the vasodilation can be assumed to be normal to the artery and, therefore, it can be modeled by only a scaling factor in that direction. Then, the vasodilation model

236		Frangi, Laclaustra, and Yang
	Table 5.1: Different similarity measures
	SSD	Sum of squared differences
	CC	Cross correlation
	GCC	Gradient image cross correlation
	JE	Joint entropy
	MI	Mutual information
	NMI	Normalized mutual information

is a similarity transformation with four degrees of freedom:
x	=	cos θ	−sy · sin θ				x		tx	(5.3)
y		sin θ	sy	·	cos θ		· y		+ ty

5.3.2.2 Registration Measure

Several registration measures have been traditionally used in medical image matching and they main ones are listed in Table 5.1.

Among the different similarity measures, normalized mutual information (NMI) is selected in this work because of its low sensitivity to the size of the overlap region [21] and higher accuracy (see section 5.4.2.2). This measure, is defined as

NMI( A, B)	=	H( A) + H(B)			(5.4)
			H( A, B)
where H( A) is the entropy of image A defined as
H( A) = − i			pi log pi		(5.5)
and H( A, B) is the joint entropy between images A and B defined as
H( A, B) = − i, j			pi, j log pi, j		(5.6)

The entropies are computed from the image histograms where pi is an approximation of the probability of occurrence of intensity value i. Similarly, the joint entropy is computed from the joint histogram where pi, j is the approximation of the probability of the occurrence of corresponding intensity pairs (i, j). Linear interpolation is used to obtain intensities in noninteger pixel values to build the joint histograms.

Computerized Analysis and Vasodilation Parameterization of FMD				237
	Table 5.2: Summary of registration parameters
	Parameter	Symbol	Value
	Registration measure	M	See Table 5.1
	No. of bins in joint histogram discretization	b	64
	Gaussian kernel width for preblurring	σ b	1
	Resampling ratio	ρ	1.5
	No. of resolution levels	r	3
	Interpolation scheme	—	Bilinear

5.3.2.3 Optimization Algorithm

A multiresolution framework proposed by Studholme et al. [22] is employed to recover the optimal transformation. The image is iteratively subsampled by a factor of two to build a multiresolution image pyramid. The registration problem is solved at each pyramid level in a coarse to fine fashion. The registration parameters found at each level are used as starting estimates for the parameters at the next level.

5.3.2.4 Summary of Registration Parameters

In Table 5.2, a summary of the parameters of the registration algorithm is provided. To compute the several registration metrics of Table 5.1, the joint histogram is discretized using 64 × 64 bins. Prior to image registration, the images are prefiltered with a Gaussian kernel of σb = 1 pixel and the images are resampled to a new pixel size of 1.5 with respect to the original size. These two steps can help to reduce small-scale noise and to reduce the computational load, and yield seemly registration results. Finally, the optimization strategy proceeds in three resolution levels. Image interpolation is carried out using a bilinear interpolation scheme.

5.3.3 Temporal Continuity

In order to perform motion compensation and vasodilation assessment, it is convenient to introduce prior knowledge about the smooth nature of the arterial vasodilation process. This a priori information could be used to filter out sharp transitions in the vasodilation parameter, which arise as a consequence

238	Frangi, Laclaustra, and Yang

of registration errors and which are not physiologically plausible. Moreover, these registration errors could easily propagate to the following frames thus invalidating all subsequent measurements.

To avoid error propagation and impose constraints on the vasodilation dynamics, a recursive filter is employed in order to improve the estimation of the initial registration parameters. In the next two subsections, the elaboration of the starting estimates in the motion compensation and in the vasodilation stages are presented.

5.3.3.1 Starting Estimate in Motion Compensation

The motion compensation phase involves three parameters: two translations, tx and ty, and a rotation angle, θ . A parameter vector is thus defined by x = {tx, ty, θ }. Because of the strong nonstationary behavior of motion artifacts, it is not possible to derive an elaborated linear model of the dynamics the parameter vector. Therefore, a simple first-order auto-regressive model, AR(1), was assumed to predict a suitable initialization for the registration algorithm in the nth frame, xˆ (n). This was done according to

xˆ (n) = xˆ (n − 1) + γ (x(n − 1) − xˆ (n − 1))

(5.7)

where x(n) denotes the parameter vector output after registering the nth frame. Note that there is some implicit delay between xˆ (n) and x(n) since they refer to parameter vectors before and after the registration process. Equation (5.7) introduces a systematic inertia to changes in the parameter values through the constant γ . This filtering tries to avoid falling into local minima in the parameter space during registration, which would not be temporally consistent with previous history of arterial motion. On the other hand, it might also slow down the ability to track sudden transitions coming from true motion artifacts. A value of γ = 0.1 has been empirically shown to be a good compromise between these two competing goals and it was used throughout our experiments.

5.3.3.2 Starting Estimate During Vasodilation Assessment

In this stage we assume that the translation and rotation parameters were correctly recovered at the motion compensation stage. Therefore, only the scale factor, sy, will be tracked over time using a simplified Kalman filtering scheme [23].

Computerized Analysis and Vasodilation Parameterization of FMD

239

Figure 5.5: Kalman gain. (a) Estimated σw (n) used for the computation of K (n) during vasodilation assessment. It contains the expected vasodilation dynamics along the sequence. Instants with higher value of σw (n) correspond to higher uncertainty about the chosen dynamic model and, consequently, where it has to be relaxed to accommodate for possibly sudden transitions. (b) Corresponding Kalman gain for three different measurement noise, σv (n), values corresponding to the minimum, median, and maximum noise levels, respectively, in our sequence database.

Let us assume that sy(n) can be also modeled as an AR(1) process

sy(n) = α sy(n − 1) + w(n)

(5.8)

where w(n) is white noise with variance σw , and 0 < α < 1 is the coefficient of the AR(1) model, which was chosen as α = 0.95. The scaling factor has a nonstationary behavior and, therefore, σw (n) actually changes widely over time (cf. Fig. 5.5).

Let the measurement model be

s˜y(n) = sy(n) + v(n)

(5.9)

where s˜y(n) are noisy measurements of vasodilation at frame n (obtained via image registration), and v(n) is white noise, uncorrelated with w(n). Under the assumptions of this model, it can be shown [24] that the Kalman filter state estimation equation to predict the vasodilation initialization for the registration of the nth frame, sˆy(n), is

sˆy(n) = α sˆy(n − 1) + K (n) [s˜y(n − 1) − α sˆy(n − 1)]

(5.10)

where K (n) is the Kalman filter gain. Owing to the standardized acquisition protocol, the vasodilation time series has a characteristic temporal evolution,

240	Frangi, Laclaustra, and Yang

which can be exploited to give an a priori estimation of σw (n). The value of σw (n) is high when vasodilation is expected and, it is low when no variations in the artery diameter should be found (e.g. at baselines). On the other hand, the observation noise power σv (n) will be considered constant, and it is estimated from the first 60 sec when vasodilation is known to be zero. The measurement noise is assumed to be stationary as it mainly depends on the image quality, which can be considered uniform over time for a given sequence.

The temporal evolution of σw (n) is shown in Fig. 5.5. It has been estimated from the analysis of 50 vasodilation curves (from the dataset described in section 5.2) that were free from artifacts and obtained with the computerized method but using a fixed K (n) = 0.1 to assess the vasodilation. The plot indicates the average instantaneous power over the 50 realizations. The 60 initial frames are processed with K (n) = 0 to calculate σv in each sequence, because no vasodilation is expected during this interval and any variation here can be regarded as measurement noise.

The values of σw (n) and σv determine the Kalman gain, K (n), using the following equation [24]:

K (n)	=		α2 K (n − 1) + σw2 (n)/σv2							(5.11)
		α2 K (n		−	1)	+	w	v +	1
							σ 2	(n)/σ 2

5.4 Computerized FMD Analysis Performance

5.4.1 Examples

Image registration between two frames searches for the transformation that puts them into correspondence. To visually illustrate the algorithm performance, four examples are shown in Fig. 5.6. In four sequences, the reference frame has been aligned with the frame showing maximum FMD.

5.4.2 Evaluation

Three properties of the proposed method are analyzed: accuracy (agreement with the gold standard), reproducibility (repeatability), and robustness (degree of automation of the measurement without apparent failure). To evaluate

Computerized Analysis and Vasodilation Parameterization of FMD

241

Figure 5.6: Four registration examples of four sequences of different image quality. Each figure shows the flow-mediated dilation frame where maximum vasodilation occurs (left half) registered to the reference B1 frame (right half).

accuracy, a number of expert manual measurements are analyzed and utilized as gold-standard measurements. As well, to get a pattern to compare our method with, accuracy and reproducibility of these manual measurements are also calculated.

We define %FMD as measurement unit for all the sequence values analyzed and represent the relative diameter of a given frame in the sequence to the mean diameter over phase B1 expressed as percentage.

Two statistical methods are used. Firstly, we used Bland–Altman plots [25], a classical method to define limits of agreement between two measurement techniques as indicated by d ± 1.96SD where d is the mean difference (bias) and SD is the standard deviation of the differences.

Secondly, we used analysis of variance to estimate the variability (reproducibility) of repeated measurements on every frame. We expressed these also as coefficient of variation (CV), obtained from the mean value (m%FMD) and the standard deviation (SD%FMD) of the %FMD measurements as indicated

CV =	SD%FMD	(5.12)
	m%FMD