Ultrasonic mapping of the microvasculature: signal alignment

Ultrasound in Med. & Biol., Vol. 24, No. 6, pp. 809 – 824, 1998 Copyright © 1998 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/98 $19.00 1 .00

PII S0301-5629(98)00042-8

● Original Contribution ULTRASONIC MAPPING OF THE MICROVASCULATURE: SIGNAL ALIGNMENT B. G. ZAGAR,† R. J. FORNARIS,‡

AND

K. W. FERRARA‡*

†

Department of Electrical Engineering, Technical University Graz, Graz, Austria, and ‡Department of Biomedical Engineering, University of Virginia, Charlottesville, VA (Received 19 September 1997; in final form 10 March 1998)

Abstract—The ultimate goal of this work was the development of a system capable of estimating the low flow velocities in the microvasculature. Estimation of low velocity flow within these vessels is challenging due to the small signal levels and the effect of cardiac and respiratory motion. Realignment of the signal from a single line-of–sight to remove physiological tissue motion is a critical part of the process of small-vessel flow mapping, and our methods for this alignment are considered in this paper. Each method involves the correlation of pulses acquired from the same line-of–sight. The first method involves the correlation of adjacent pulses (nearestneighbor), the second involves a single reference line and the third involves averaging the correlation over a set of reference lines. We find that a nearest-neighbor strategy is suboptimal, and that strategies involving a global reference line are superior. A bound on the variance of estimates of the location of the correlation peak is presented. This bound allows us to consider our results in comparison with an absolute limit. Finally, a new algorithm allowing for alignment between lines-of–sight is described, and initial results are presented. Such an algorithm does, in fact, reduce jitter, correct for tissue motion and enables us to better visualize vessel continuity. We find that vessels as small as 40 mm can be mapped in two dimensions using a 50-MHz transducer. © 1998 World Federation for Ultrasound in Medicine & Biology. Key Words: Signal alignment, Cross-correlation, Velocity estimation, Blood flow mapping.

INTRODUCTION

signal alignment to remove motion induced by cardiac, respiratory and other involuntary sources. First, physiological tissue motion velocities have been found to resemble closely those of slow blood flow, necessitating motion correction before blood velocity estimation. Although the velocity of slowly moving blood may be estimated without correction for physiological motion, in practice this has not been successfully accomplished. Therefore, because these sources of movement may have the same velocity as blood flow in the microcirculation, additional parameters must be incorporated into the analysis to differentiate between them. Specifically, we used the differing spatial extents of motion due to blood flow and physiological processes; the latter typically demonstrating a uniform rate over a large depth. Thus, alignment using pulse-to–pulse correlation over a wide spatial window can remove tissue motion without disturbing the echoes from blood. In addition, echoes from tissue have higher amplitudes than those from blood. Therefore, correlation techniques that use a significant depth window and heavily weigh the signal from tissue might successfully estimate and correct for gross tissue motion.

Ultrasound provides a noninvasive technique to visualize and study tissue structures. We have previously reported the development of a high-resolution algorithm for velocity estimation (Ferrara et al. 1991, 1996b) and have now extended this work to map blood flow in small vessels in vivo. No alternative methods currently exist for mapping the velocity of blood flow in the microvasculature of opaque tissues, which is particularly challenging because cardiac and respiratory motion can confuse the analysis. Using ultrasound with center frequencies up to 50 MHz, reproducible velocity estimates as low as 0.2 mm/s, combined with an effective spatial resolution on the order of 40 mm, can be obtained (Lockwood et al. 1991). This very high resolution, both in space and velocity, increases the system’s motion sensitivity, ultimately affecting blood velocity estimation. In this paper, we concentrate on the problem of

*Address correspondence to: Katherine W. Ferrara, Department of Biomedical Engineering, Box 377, Health Science Center, Charlottesville, VA 22908, USA. 809

810

Ultrasound in Medicine and Biology

To provide an absolute basis for alignment performance comparison, an expression for the variance of the location of the correlation peak is first derived. Evaluating this expression for the experimental conditions, we found that it predicts a residual jitter after alignment on the order of 1 ns. Then several correlation-based alignment methods are evaluated, and we demonstrate that they result in a considerable subjective improvement in pulse-to–pulse alignment. Using data from a specially constructed phantom and two different measures of residual jitter, we have found that the use of a single reference line is superior to a nearest-neighbor algorithm. In addition, we show that performing the alignment following a data interpolation algorithm producing 1 ns pixels is nearly optimal. One long-term objective of this research is to combat major causes of blindness by providing quantitative and comprehensive assays of perfusion in ocular structures. The unique information provided by such noninvasive techniques may be crucial for diagnosing ocular diseases, assessing their severity, and monitoring treatment efficacy. The techniques could also provide key information regarding normal perfusion levels, as well as physiological and age-dependent variations. The importance of such data has long been acknowledged but has previously been unattainable due to the small sizes of ocular components, the presence of overlying opaque tissues, and the need to employ noninvasive techniques within delicate ocular structures. Of particular significance are studies of glaucoma, a major cause of blindness whose incidence is expected to rise in the USA as the population ages. Our high-resolution techniques allow us to characterize ciliary body perfusion in specific subcategories of glaucoma and monitor changes induced by pharmacological agents, surgery and laser therapy. The techniques could also characterize the course of changes in the peripheral retinal vasculature arising from elevated intraocular pressure. Experimental studies and clinical observations indicate that increased intraocular pressure alters structural and neural elements of the optic nerve head. The increased pressure does not explain other observations in chronic primary open-angle glaucoma and in normaltension glaucoma. Changes in optic nerve perfusion do occur in these cases, and a noninvasive estimate of perfusion may be of particular value. Vascular changes have also been documented in malignant melanoma. Histopathologic studies of the microcirculation architecture of uveal melanoma by Folberg et al. (1992, 1993), and Rummett et al. (1994) have demonstrated a correlation between vascular features and tumor lethality. Less et al. (1991) have shown that the microvascular architecture of solid mammary tumors differs distinctly from normal tissues. A vascular net-

Volume 24, Number 6, 1998

work exists that exhibits fluctuations in both the diameter and length of the vessel with increasing branch order. Using ultrasound parameter imaging to produce maps of effective scattering element size in uveal melanoma, Coleman et al. (1990, 1991) have shown a spatial association between vascular network patterns and intratumor regions of smaller scattering elements. We have begun an evaluation of the vasculature of human prostate tumors hosted by mouse models. With the development of engineered mouse tumor models in which human tumors grow in a superficial region of a mouse, the use of high-frequency ultrasound may be able to provide new information as to the location and size of vessels and the velocity of blood flow within the mass. Tools for mapping local vessel density throughout an intact tumor are not currently available; therefore, local vessel counts obtained through histology are used for comparison. Folkman et al. (1971) recognized the importance of tumor vascularity and hypothesized that the increased cell population required for the growth of a malignant tumor must be preceded by the production of new vessels. Recently, the mean vascular density has been shown by Bigler et al. (1993) to be significantly higher in cancerous prostate tumors with a mean vessel count of 135.8 vessels/mm2 and standard deviation of 49.1 vessels/mm2. In benign prostate tumors, their corresponding vessel count indicates a mean of 71.9 vessels/ mm2 and standard deviation of 30.9 vessels/mm2. Measurements by Baish et al. (1996) indicate that the vascular patterns of tumors are irregular, with regions of increased and decreased density having a spatial extent on the order of 600 mm. The most effective tool for the assessment of microvascular perfusion in optically accessible tissue has been laser Doppler, and these studies have shown the expected decrease in flow rate predicted from the tortuous architectures. Foltz et al. (1995) have shown that laser Doppler flow measurements are significantly decreased in tumor microvasculature in comparison to control sites. Other evidence of this decrease includes simulations of flow based on known vessel architecture and resistance by Netti et al. (1996), which predict lower flow rates and the potential for stasis in the tumor center. Measurements of vascular resistance in excised tumors from previous studies (Less et al. 1998) show a significant increase in resistance in tumors (83 in comparison with muscle) and, therefore, also predict a decrease in microvascular flow rate. First, an overview of the signal processing required for mapping low-velocity blood flow is presented. The experimental system used to acquire this data is then summarized and sample in vitro and in vivo data sets are presented. An expression for the variance in the location of the correlation peak is then derived. Three specific

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

algorithms for alignment within a line-of–sight are detailed and their resulting performances evaluated, and an algorithm for alignment between lines-of–sight is also considered. Finally, the resulting blood-flow map following alignment is presented.

811

SIGNAL-PROCESSING METHODS

tude of the signal from the vessel wall was evaluated and shown to be 12–18 dB above the signal from blood at the center frequency of 38 MHz (Ferrara et al. 1996b). A 4th-order Butterworth IIR step-initialized wall filter producing more than 18 dB of attenuation was used in these feasibility studies, producing reproducible velocity estimates as low as 0.2 mm/s.

The critical obstacle in mapping the microvasculature is differentiating slowly moving tissue from blood flow. The signal-processing methods employed were optimized to produce high spatial and velocity resolutions. Our signal processing involves the following steps: alignment within and between lines-of–sight, the application of initialized IIR wall filters, velocity estimation that involves tracking small groups of red blood cells, the application of a threshold algorithm to differentiate flow from noise and post-processing of the velocity profile. Alignment methods are described in detail later, and other signal-processing algorithms are briefly summarized below. Two measures of alignment performance are evaluated, the first being the magnitude of the pulseto–pulse correlation and the second being the postwallfiltered power. In a region of stationary tissue without resolvable blood vessels, the pulse-to–pulse correlation should be maximized by proper signal alignment. Also, in such a region, the signal that remains after the application of a high-pass filter (the postwall-filtered power) is typically produced by high-frequency noise or misalignment due to jitter or cardiac and respiratory motion. Thus, in such data sets, we evaluate the postwall-filtered power for a typical wall filter and use this as an indication of the quality of the alignment. In a region with resolvable blood vessels, the postwall-filtered power is significant and can be compared to a predetermined threshold to reject noise.

Velocity estimation In previous studies (Ferrara and Algazi 1991a,b; Ferrara et al. 1996b) we presented a wideband point maximum likelihood estimator (WMLE), which uses both the change in position of the tracked scatterers and the Doppler shift of the returned echo with respect to the excitation. Both are produced by the movement of small groups of red blood cells. Results demonstrate that the WMLE allows for significantly faster velocity estimation, as compared to a cross-correlation estimator, and yields a higher signal-to–noise ratio with a specific model. The particular signal model used here assumes the presence of a single dominant velocity component at each location at a given time. The WMLE for a discrete velocity value, y, corresponding to a two way travel time, t (d), is given by the maximum over y of the likelihood, l( y ), as shown in Equation (1). The presence of y within the envelope, s90 ( . . . ), and the exponential matches the change in delay and frequency shift caused by the movement of the red blood cells, respectively. Derived using the slowly fluctuating point target approximation, the likelihood involves the sum over the pulse train, indexed by k, and the integral over the temporal axis, t. The integrand is the product of the complex envelope of the received signal, r9(t); an estimate of the complex envelope of the received signal, denoted by s90 ( . . . ); and the Doppler shift. The WMLE is given by:

Wall filter A wall filter is a high-pass filter that removes the signal from stationary or slowly moving tissue to allow for accurate blood velocity estimation. Both step-initialized and exponentially-initialized IIR wall filters (Kadi and Loupas 1995; Ferrara et al. 1996b) effectively suppress by 18 dB or more any signal components caused by scatterers with velocities below '0.1 mm/s. However, these wall filters require the correction of gross tissue motion and pulse-to–pulse jitter before filtering, because even mild misalignment causes the filter to ring. The velocity estimator then interprets this ringing as blood flow and generates erroneous velocity estimates, altering resulting flow maps and preventing accurate 3-D reconstruction. Step initialization was chosen for this application due to the low velocity of this wall motion. The ampli-

OU O

P221

1 l~y! 5 P2

l50

E

1`

3

2`

1 P1

~l11! P121

exp@ j2vcykT/c#

k5lP1

U

2

r9~t!s90~t 2 t~d! 2 kT@1 1 2y/c#!dt ,

(1)

where v c indicates the angular frequency, c indicates the acoustic propagation velocity, T indicates the pulse repetition period, P 1 indicates the number of pulses combined coherently, and P 2 indicates the number of functions combined incoherently to form the likelihood. This strategy significantly decreases the variance of the estimate, such that it remains constant over a wide range of velocities. The WMLE tracks the spatial volume in time and

812

Ultrasound in Medicine and Biology

Volume 24, Number 6, 1998

Table 1. Parameters characterizing the experimental setup Parameter

Human

Mouse

Rabbit

Phantom

26 dB bandwidth B in MHz Transducer center frequency in MHz Center frequency of echoes f 0 in MHz Sampling frequency in MHz Pulse repetition frequency in Hz Number of pulses per line-of–sight

22 50 37 250 1000 16

12 35 27 500 200 30

22 50 37 250 1000 30

18 50 25 250 1000 100

space, and the observation time of a group of scatterers is not limited by the signal length. Therefore, a very short duration signal can be transmitted to improve spatial resolution. Threshold and postprocessing After velocity estimation at each axial location using the WMLE, thresholds for the magnitude of the likelihood and postwall-filtered power are used to differentiate between velocity estimates produced by noise sources and estimates corresponding to regions of blood flow. These estimates are then postprocessed using an adaptive nonlinear filter and a measure of 3-D continuity of detected blood vessels to further reduce noise and improve accuracy. Additional details on these algorithms have been presented by Ferrara et al. (1996a).

beam width is 35 mm. A 35-MHz center frequency transducer used for the tumor study is the PI35-2-F1.00. A LeCroy (Chestnut Ridge, NY) digitizing oscilloscope was used to acquire the signal and apply a gate to restrict the depth of interest. One thousand samples were acquired per pulse at a sampling rate of 250 MHz. In each case, the high-frequency ultrasound transducer was coupled to the tissue under evaluation through a water bath, and a 3-D region evaluated by mechanically scanning the transducer over the region. An independent calibration of the 50-MHz transducer/pulser combination determined an ISPTA intensity of 0.65 mW/cm2 for a transmitted energy of 2 J and pulse-repetition frequency of 1000 Hz, parameters used in our study. SAMPLE DATA SETS

EXPERIMENTAL SYSTEM We report on two sets of in vivo experiments in this paper involving the eye and a tumor model. The parameters used in these studies are summarized in Table 1. A more detailed description of our experimental system has been previously presented (Ferrara et al. 1996b). For the ophthalmic study, the higher center frequency of 50 MHz was chosen, together with a slightly higher pulse repetition rate, due to the higher velocities of flow. Fewer pulses per line-of–sight were transmitted in the evaluation of the human eye because higher velocities are detected and can be discriminated from tissue motion more easily. In the evaluation of the mouse tumor models, a lower center frequency of 35 MHz was chosen to improve the depth of penetration, and a pulse repetition frequency of 200 Hz was chosen due to the very low flow velocity expected within small vessels in tumors. The experimental system consists of a Panametrics 5900 pulser/receiver with a 200-MHz bandwidth providing impulse excitation of a PVDF transducer. The received signal is amplified by 54 dB within the pulser/receiver. The 50-MHz center frequency transducer used is the Panametrics (Waltham, MA) PI50-2-F0.50, with a 15-dB transducer bandwidth from 10 –50 MHz. The measured center frequency in the returned echoes with this transducer is 38 MHz, the one-sided 6 dB lateral

Four sample data sets are introduced here and are used in the following sections to demonstrate the magnitude of the distortion that can be introduced by physiological motion. First, we consider an M-mode image of a rabbit eye, obtained from an anesthetized rabbit. Significant blood flow is demonstrated, and this sample image will be used to demonstrate the performance of the alignment algorithms. Second, a scan of the human eye is presented to demonstrate the jitter produced when the subject attempts to fix his vision on a light source. Third, a set of phantom data not having blood flow was constructed and will allow for quantitative evaluation of algorithm performance. Finally, we present data from a mouse model in which no resolvable blood vessels are present. The severe motion in this image can be removed and this case will again be used quantitatively to compare the results of signal alignment. For all in vivo data sets, each line-of–sight is represented by a single M-mode image, with a stacked set of echoes from individual pulses along the vertical axis. In each case, sets of black and white stripes have been placed along the vertical edges of the image to indicate the lines-of–sight. All pulses within a single line-of–sight are shown within this stripe. Additional M-mode images from subsequent lines-of–sight are then stacked to form a 2-D set of data. In each M-mode image, stationary

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

813

Fig. 2. An enlarged subimage from Figure 1 demonstrating tissue motion in the ciliary body from a rabbit eye. Fig. 1. M-mode image of the anterior segment of a rabbit eye.

targets appear as vertical stripes; moving targets appear as slanted lines whose slope indicates the magnitude of the velocity. Figure 1 shows the M-mode image of the anterior segment of a rabbit eye, scanned with the PI50-2-F0.50 transducer and sampled at 250 MHz, with a pulse repetition frequency of 1 kHz. Forty lines-of–sight, spaced 40 mm apart, with 30 pulses per line-of–sight, were acquired. Three large regions of flow are visible, with one in the upper region of the image and two in the lower region, each surrounded by an ellipse. In the upper region, flow within a large vessel is visible at the central depth bins. In the lower region, we observe two closely spaced blood vessels with diameters of 200 mm and 40 mm and with blood flow in opposite directions. The two vessels are 75 mm apart. We use this image to demonstrate the effects of the alignment algorithms and the resulting blood-flow maps. In the absence of cardiac and respiratory motion, echoes from stationary tissue would be vertical. In this case, they are slightly slanted or curved. A small segment of this image is also shown in Figure 2, where pulse-to–pulse jitter is more obvious. Figure 3a shows an M-mode image of the anterior segment of a human eye, acquired with the PI50-2-F0.50

transducer. Shown here are 32 lines-of–sight 40 mm apart; 16 pulses per line-of–sight were acquired, with a pulse-repetition frequency of 1 kHz. The rapid contractions of the ciliary muscles further complicate alignment by introducing extensive tissue motion. This motion produces highly slanted echoes, as well as decorrelation due to out-of–plane motion. Figure 3b shows a single M-mode line-of–sight from a sponge phantom. In this case, the transducer was rapidly moved toward and away from the sponge at a variable velocity to simulate the effects of respiratory and cardiac motion. The instantaneous velocity is approximately 1 mm/s, similar to that of the in vivo images. Figure 3c shows a portion of an M-mode scan of a poorly vascularized region in a mouse tumor model, insonified with a 35-MHz transducer. For each line-of– sight, 30 pulses, with 500 samples per pulse, were acquired using a sampling frequency of 500 MHz and a pulse repetition frequency of 200 Hz. All three in vivo images show both low-velocity gross tissue motion and less visible, but still severe, pulse-to–pulse jitter, both of which obscure blood flow. The mean velocity associated with the subject’s motion, ^uyu&, and the instantaneous velocity, yinstant, were determined through spectral analysis of the necessary shift,

814

Ultrasound in Medicine and Biology

Volume 24, Number 6, 1998

cannot reduce this source of estimator variance, this knowledge sets the lower bound on residual correlation error. Also, this timing jitter advises a limit on the sampling increment in the digital interpolation used prior to these algorithms, because further interpolation is not expected to improve performance. By evaluating the performance of the alignment algorithms as a function of digital interpolation, this limit on interpolation is supported. We describe each of the three algorithms for alignment within single lines-of–sight, as well as an algorithm for alignment between lines-of–sight. The resulting alignment of the specified data sets is presented and evaluated to measure and compare the individual performance of each algorithm.

Fig. 3. (a) M-mode image of the ciliary body of a human eye; (b) phantom data used in the comparison of alignment quality; (c) portion of an M-mode scan of a tumor mouse model.

S( p i ), to align pulses (Table 2). Although visual inspection of the M-mode images in Figs. 1 to 3 indicates that the tumor mouse model displays the highest relative shift in depth over the set of pulses acquired, it does not exhibit the fastest motion because the pulse repetition frequency is only one fifth of that used in the human and rabbit experiments. ALIGNMENT As Figures 1 to 3b indicate, gross tissue motion and pulse-to–pulse jitter can become quite severe, so the development of alignment algorithms capable of motion correction and jitter reduction is essential. The variance in estimates of the location of the correlation peak produced by noise will first be derived. Because alignment

Table 2. Mean and instantaneous velocity of tissue or phantom with respect to the transducer Parameter

Human

Mouse

Rabbit

Phantom

^uyu& in mm/s yinstant in mm/s

0.45 1.35

0.25 0.46

0.17 0.87

0.2831 2.580

Theoretical confidence interval for the location of peak correlation values Before discussing the different alignment strategies, measures of alignment performance must be developed, and few authors have reported work specifically in that field (Walker and Trahey 1995; Jensen 1996). To provide an absolute measure for performance comparison, a confidence interval for the bound on the estimate of the location of the correlation peak was first derived. Using this measure, the expected noise-induced residual jitter can be determined and the performance of our algorithms evaluated. The jitter and motion of pulse p i (t) that must be corrected is modeled as a random delay with respect to the precise occurrence of the pulse p 0 (t) delayed by iT (Figure 4). This delay occurs in addition to the round-trip delay, t (d), associated with scatterers in depth, d. The delay is drawn from a random distribution and modeled as a random variable that follows a Gaussian distribution with zero mean and variance, s2(t). In addition, each pulse, p i (t); 0 # i , P max, has an additive noise source, n i (t), associated with it. All n i (t) are assumed to be jointly Gaussian white noise sources that are uncorrelated with respect to one another and to the excitation pulse, p 0 (t 2 iT), (Equation (2)): p i~t! 5 p 0@t 2 iT 2 t i 2 t ~d!# 1 n i~t!.

(2)

Following a derivation by Bendat and Piersol (1986) and assuming ergodicity of the model (see Appendix), the location of the peak in the cross-correlation function, R p r,p i( t ), for pulse p i (t) and the reference pulse, p r (t), exactly determines relative delays, ( t i 2 t r 1 (i 2 r)T), between pulses that must be corrected during the alignment procedure (Equation (3)): R pr,pi~ t ! 5 R p0,p0~ t 2 @ t i 2 t r#!.

(3)

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

815

Fig. 4. Model describing the introduction of pulse-to–pulse jitter and motion.

Because each pulse is limited in duration ( p 0 (t) Þ 0 for 0 # t , t max), only an estimate of the autocorrelation function, Rˆ p 0,p 0( t ), can be determined from acquired data (Equation (4)):

Rˆ p0,p0~ t ! 5

1 t max 2 t

E

Given the above assumptions, the normalized RMSerror «{Rˆ p 0,p 0(0)} can be determined to be (Bendat and Piersol 1986):

tmax2t

p 0~t! z p 0~t 1 t !dt,

var$Rˆ p0,p0~0!% < «$Rˆ p0,p0~0!% 5 R p20,p0~0!

0

2 T

E

1`

R p20,p0~ j !d j

2`

R p20,p0~0!

. (6)

0 # t , t max.

(4)

Because the alignment algorithms discussed later rely upon determining the exact location of the peak delay, ( t i 2 t r ), of the cross-correlation function, we must derive an expression for the variance of that random variable to quantify alignment quality. This also sets a reasonable limit to the amount of interpolation applied before velocity estimation is performed. Assuming that the spectral content of the pulse, p 0 (t), is band-pass white noise with bandwidth, B, centered at frequency, f 0 , and that the uncorrelated noise sources, n i (t); 0 # i , P max, are low-pass white noise sources, the variance, s 2 ( t i 2 t r ), of the location of the peak in the correlation function can be determined as (see Appendix for derivation):

s2~ti 2 tr! 5

6/Î3 z «$Rˆp0,p0~0!%, ~pB!2 1 3~2pf0!2

(5)

In further discussions on the error bounds for alignment algorithms, the normalized RMS-error, «{Rˆ p 0,p 0(0)}, will be directly estimated from the recorded data. The four M-mode scans shown in Figs. 1 to 3b were statistically analyzed to calculate the standard error «{Rˆ p 0,p 0(0)} directly from a set of '500 pulses using the middle term in Equation (6). This error figure was then used to calculate s 2 ( t i 2 t r ) in Equation (5), shown in Table 3. This leads to the conclusion that, for all sample data sets analyzed in this paper, a reasonable limit in the sampling rate 1/T is on the order of 1 ns independent of whether this sampling rate is achieved by digital interpolation or used within data acquisition. Algorithms for alignment of an array of pulses acquired along a single line-of–sight The variance in estimates of the location of the correlation peak was derived in the previous section for

with B spectral bandwidth of pulse p 0 (t) f0 center frequency «{ . . . } normalized root mean square (RMS)-error of {...} Rˆp0,p0(0) autocorrelation estimate of pulse p 0 (t) at lag t50

Table 3. Predicted jitter in the correlation peak for different subjects Parameter

Human

Mouse

Rabbit

Phantom

«{Rˆ p 0,p 0(0)} s ( t i 2 t r ) in ns

0.0326 0.822

0.0327 1.13

0.0152 0.561

0.022 0.993

816

Ultrasound in Medicine and Biology

continuous time data. In this section, discrete time algorithms are detailed for the case in which a small array of radiofrequency signals (pulses) is obtained from a single line-of–sight. The result of all discussed algorithms is the vector, S( p i ), indicating the shift necessary to align pulse p i with the other pulses within a single line-of– sight. The shift is given in integer multiples of the sampling increment, which has been digitally interpolated in most cases. All algorithms assume that stationary targets, such as vessel walls, should not move with respect to the transducer and, therefore, their returns should be perfectly aligned, yielding a pattern of vertical structures. These algorithms will vary in their response to decorrelation of the signal within a line-of–sight. Each of these algorithms is based on the assumption that the location of the cross-correlation peak of two pulses containing a common signal (a delayed version of pulse p 0 (t)) will give an unbiased estimate of the delay between these two pulses. The cross-correlation estimate Rˆ p i,p r[m] for discrete time data is given in Equation (7) (Marple 1987). Based upon this cross-correlation approach, various pulse combinations are used to better align the signal returns. Motion of the target across the transducer beam would reduce the similarity between pairs of pulses and, therefore, lead to decorrelation.

5

O

1 p ~n 1 m! pr~n!, N 2 umu n50 i 0#m,N21 Rˆpi,pr@m# 5 1 N2umu21 p ~n! pr~n 1 umu!, N 2 umu n50 i 2 ~N 2 1! , m , 0 N2m21

O

(7)

with: Rˆpi,pr[m] m pi(n)

Volume 24, Number 6, 1998

Fig. 5. The central portion of the cross-correlation function Rˆ p i,p r[m], as defined in Equation (7) for pulses 60 to 80 of the data shown in Figure 3b. Lag shown in 4 ns increments.

within a line-of–sight as a reference pulse and calculates the relative shift of all other pulses with respect to that reference pulse.

H

J

S I~ p i!u pr 5 arg max $Rˆ pi,pr@m#% . m

(8)

S I ( p i ) specifies the relative shift of pulse p i with respect to the reference pulse, p r , which is kept constant within a line-of–sight. This strategy, being the simplest of the three to be described, works well if no significant decorrelation from pulse to pulse alters the internal structure of the ultrasound return. It also performs faster than the others. If, however, stationary scatterers move out of the focal volume of the transducer due to the subject’s motion, significant decorrelation over the set of pulses can

normalized cross-correlation estimate of the returns for pulses i and r time lag argument of the cross-correlation estimate ultrasound return of pulse p i ; 0 # i , P max with index n denoting depth

Figure 5 shows the central part of a particular crosscorrelation function estimate, Rˆ p i,p r[m]. The pulses to be aligned were taken from the phantom data set shown in Figure 3b (pulses 60 to 80). The crest of this function is indicated by a thick yellow line, and its location represents the relative shift of pulse p i with respect to the reference pulse, p r , pulse number 60 in Figure 3b. Appropriately applying the negative of the determined shift to each of the pulses yields a nearly perfectly aligned image. Single reference line algorithm. The first strategy, a reference line algorithm, selects a single pulse, p r (n),

Fig. 6. Aligned data from a small region of phantom in Figure 3b.

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

occur, thereby impairing the quality of alignment (see Figure 6, second image). Incremental algorithm. The second strategy, an incremental or nearest-neighbor algorithm, is able to cope with slowly decorrelating returns from consecutive pulses and results in an accurate alignment. The relative displacement DS( p j , p j21 ) of consecutive pulses from ( j 2 1) 3 j is determined by finding the maximum of the cross-correlation estimate over the time-lag argument, m, as given in Equation (9):

H

J

2 # j # P max.

5 (9)

The total pulse delay, S II ( p i ); 1 # i # P max, to be applied to each of the pulses is the accumulated displacement, DS( p j ), up to index i, as defined in Equation (10):

SII~ pi! 5

HO

0,

i51

i21

j51

DS~ pj!,

2 # i # Pmax

(11), the reference line algorithm, Equation (8) is used repetitively so that each pulse within a line-of–sight serves as a reference pulse once. If quantization errors and decorrelation effects are absent, all shift vectors S I ( p i )u @p r would be equal. Because this is not the case, ensemble averaging over all reference pulses is capable of both reducing effects of decorrelation and averaging over quantization errors. S III ( p i ) is the ensemble average ^ . . . & i over all relative shifts for all pairs of pulses ( p i, p ); p i Þ p j . j S III~ p i! 5 ^S I~ p i!up j& 1#j#Pmax

DS~ p j! 5 arg max $Rˆ pj,pj21@m#% , m

817

(10)

This strategy performs far better in the presence of strong decorrelation than the reference line algorithm because the decorrelation between adjacent pulses is very limited, in most cases. However, because this algorithm compares adjacent pulses, it does not have any global perspective; thus, it cannot compensate for mild misalignment from pulse to pulse. This problem is especially dramatic because all computations are performed on discrete time data, so quantization effects can become significant. In the worst case, if the relative shift from pulse j 2 1 to pulse j for all j is less than one half of the sampling increment, the algorithm will yield all DS( p j ) to be zero, yielding a cumulative sum, S II ( p i ), up to pulse i of zero. Figure 6 shows a small region of the phantom data set, with the prealigned data in the top region of the image. An example of the alignment attainable with the incremental algorithm is shown in Figure 6 in the third image. Within the phantom data set shown in Figure 6, the motion of the transducer was nearly perpendicular to the phantom and, therefore, the data remains highly correlated over a large ensemble of pulses. A diagonal line runs through the center of the image, and was produced by a decorrelation artifact. In this case, the quality of the alignment is primarily limited by quantization effects, where this is observed by the remaining jitter. Averaging algorithm. A novel scheme, the averaging algorithm, which is an improved version of the reference line algorithm, is introduced here. According to Equation

O argH max $Rˆ

Pmax

1 P max

j51

J

@m#% ,

pi,pj

m

1 # i # P max.

(11)

The performance of the averaging algorithm is subjectively superior to the previous algorithms, and the resulting aligned data are presented in the bottom image of Figure 6. In the following section, the performance of this algorithm will be shown to be quantitatively superior to the previous algorithms, but this performance is achieved at the price of a greatly increased numerical complexity. Alignment between lines-of–sight Alignment within single lines-of–sight improves blood velocity estimation and allows for more accurate vessel reconstruction. Alignment between lines-of–sight is now considered as a method to further reduce the effects of tissue motion, assess vessel continuity and improve flow mapping and 3-D vessel reconstruction. Data acquisition along a set of discrete lines-of–sight introduces spatial discontinuities; however, slow tissue motion aggravates these discontinuities and greatly affects the resulting flow map. Using a sufficiently small spacing between lines-of–sight, resulting discontinuities are primarily attributable to tissue motion. The unaligned data from the anterior segment of a rabbit eye is shown in Figure 1. Because the stationary scatterers, represented by the black and white bands, are not vertical, jitter reduction and motion correction within and between lines-of–sight are necessary. Again, alignment between lines-of–sight is based upon the standard cross-correlation between two pulses, chosen to represent entire lines-of–sight.

Rˆli,lj@m# 5

5

1 N 2 umu

O

N2m21

li~n 1 m!lj~n!,

n50

0#m#N21

1 N 2 umu

O

N2umu21

li~n!lj~n 1 umu!,

n50

2~N 2 1! # m , 0

(12)

818

Ultrasound in Medicine and Biology

Volume 24, Number 6, 1998

Fig. 7. Averaged cross-covariance vs. distance between pulses for phantom data.

li Last pulse from line-of–sight i lj First pulse from the next line-of–sight j for j . i m Time lag argument of the cross-correlation estimate Lmax Number of lines-of–sight Rˆli,lj[m] Normalized cross-correlation estimate of the returns for the last pulse and first pulse of their respective lines-of–sight l i (n) and l j (n), where 1 # i # L max 2 1 and j . i By extending the incremental algorithm, the relative shift between the last pulse of one line-of–sight and the first pulse of the subsequent line-of–sight can be determined. The entire block of pulses within the M-line is then shifted by the determined amount, thereby aligning linesof–sight.

RESULTS To evaluate the performance of the alignment algorithms, we first examined realignment of the test images for a subjective evaluation of performance, and then considered two quantitative measures. Finally, we show a velocity map attained from the interpolated and aligned data set corresponding to Figure 1. Results of alignment We begin with the phantom data set in which a sponge moved with respect to the transducer. The cor-

Table 4. Postwall-filtered power for phantom data set (no fluid flow) Algorithm

Postwall-filtered power

Unaligned Incremental Reference line Averaging

156890 23185 19016 15553

Fig. 8. Results of aligning phantom data (Figure 3b). Interpolation factors used are indicated.

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

relation magnitude and postwall-filtered power are evaluated for this case. Similar measures are computed for the mouse data set that does not contain detectable blood flow or resolvable blood vessels. Applying each of the algorithms to these raw phantom data and enlarging small regions of the aligned data, as shown in Figure 6, further enables us to compare the performance of the three alignment strategies. Clearly, the averaging algorithm outperforms the other algorithms, indicated by the smoother appearance of the vertical bands. To better interpret the performance of these algorithms applied to the phantom data set, the normalized cross-covariance functions, C p i,p j(0), as a measure of linear dependence between pairs of pulses p i and p j , were calculated from the alignment result and averaged over all lines-of–sight available from the phantom data. As the distance between pulses increased, the correlation peak for all of the alignment algorithms decreased, indicating decorrelation. The main diagonal of the 2-D crosscovariance function was extracted from each of the correlation surfaces in order to better visualize the differences between the three algorithms, as shown in Figure 7. First, note that, for the prealigned data, the maximum correlation magnitude would drop to zero at an interpulse lag of 20 pulses and, thus, all three algorithms produce a substantial improvement in the correlation coefficient. Although the correlation coefficients for all three algorithms are similar, slight differences are visible where the averaging and reference line algorithms outperform the incremental algorithm. The oscillatory nature of the correlation coefficient stems from the periodicity of the RF signal that is clearly visible, although small on an absolute scale. Our second measure of alignment performance, the postwall-filtered power, indicates degrees of misalignment, where a lower postwall-filtered power corresponds to better alignment. The postwall-filtered power for the unaligned phantom data set and the results of each of the alignment algorithms are listed in Table 4. The large decrease in postwall-filtered power is attributed to improved alignment. Again, the averaging and reference line algorithms outperform the incremental algorithm, as indicated by their respective values. Alignment of digitally interpolated data. Using a typical sampling period, T, of 4 ns (equivalent to 250 MHz sampling frequency) for the data-acquisition system, the variance in the location of the correlation peak due to the residual quantization error alone would be s(tres) 5 T/ =12 5 1.15 ns, which is larger than the theoretical limit attainable for continuous time data (see Table 3). Thus, interpolation to a higher sampling rate may reduce overall variance in the location of the cor-

819

Fig. 9. Average postwall-filtered power in aligned images as a function of the interpolation factor used.

relation peak. To illustrate the benefits of implementing the alignment procedure on digitally interpolated data (Crochiere and Rabiner 1983), various interpolation factors were used to align the raw phantom data, shown in Figure 3b. To exaggerate even small residual misalign-

Fig. 10. (a) Small portion of the M-mode image of a tumor mouse model shown in Figure 3c; (b) portion of the M-mode image shown in Figure 3c after applying the reference line algorithm.

820

Ultrasound in Medicine and Biology

Volume 24, Number 6, 1998

Fig. 11. Averaged cross-covariance vs. distance between pulses for the mouse data presented in Figure 3c (solid line 5 averaging algorithm; dotted line 5 reference line algorithm; and dashed line 5 incremental algorithm).

ment, enlarged subimages are displayed in Figure 8. The top left subimage in Figure 8 exhibits the resulting alignment without interpolation, the middle left image applies an interpolation factor of 2, bottom left uses a factor of 4, top right uses a factor of 8, middle right uses a factor of 16, and the bottom right uses a factor of 32. Improved alignment due to interpolation is subjectively apparent up to an interpolation factor of four (equivalent to a sampling increment T of 1 ns), which corresponds well with the theoretical bounds that are on the order of '1 ns, as shown in Table 3. To obtain a more quantitative measure of alignment quality, the signal power after the application of a step-initialized fourth-order Butterworth filter was calculated and is displayed in Figure 9. Initially, a large reduction in postwall-filtered power occurs due to interpolation up to a factor of four. However, due to residual noise, only slight additional reductions in postwall-filtered power are observed if interpolation factors greater than four are used. Alignment of mouse data. The M-mode image of the tumor mouse model (see Figure 3c) exhibits the greatest degree of apparent motion of all the subjects discussed here. This is mostly due to the decrease in the pulse repetition rate as compared to all of the other acquisitions, as is shown in Table 3. A portion of the

mouse data is used to demonstrate the efficiency of the alignment strategies. Figure 10a depicts the input data to these algorithms. The white lines separating blocks of pulses within single lines-of–sight were introduced to clarify the viewer’s perception of vertical structures. Figure 10b shows the result of applying the reference line algorithm to the data set, where tissue echoes are now vertical within each line-of–sight. Again, the normalized cross-covariance functions, C p i,p j(0), were calculated and averaged over all lines-of– sight for the tumor mouse model. The main diagonal of the 2-D cross-covariance function was extracted from each of the correlation surfaces in order to better visualize the differences between the three algorithms (Figure 11). Although the correlation coefficients for all three algorithms are similar, slight differences are visible where the averaging and reference line algorithms outperform the incremental algorithm. As the distance between pulses increases, the incremental algorithm’s correlation coefficient decreases more rapidly than those of the other algorithms, which perform similarly. Comparing Figures 7 and 11, we find that white noise causes a rapid decrease between lag 0 and 1 in the correlation magnitude for in vivo data. This is most likely due to a higher amplifier gain required for the in vivo

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

821

Table 5. Postwall-filtered power for mouse tumor model without detectable blood flow Algorithm

Postwall-filtered power

Unaligned Incremental Reference line Averaging

146200 6705 5814 5530

experiment. Second, we note that the in vivo correlation does not rise at large lags because there is a greater spatial variation in the scattering structures of tissue. Finally, we note that the decorrelation is slightly greater for the in vivo data, most likely due to noise and lateral physiological motion. The second quantitative measure of alignment algorithm performance, the postwall-filtered power, was then determined for the mouse data set. Again, since the postwall-filtered power indicates the degree of residual misalignment, a decrease in power corresponds to better alignment, and the application of our algorithms does, in fact, reduce the postwall-filtered power, which is presented in Table 5. Although all three algorithms drastically reduce the postwall-filtered power, the averaging and reference line algorithms outperform the incremental algorithm, for this case as well as our phantom data results. Alignment of structures between lines-of–sight. In Figure 12, the same ciliary body shown in Figure 1 has been aligned both within and between lines of sight. Two regions of flow, faint and obscured by jitter and tissue motion in the first image, are more easily distinguished. Also, even though we have aligned between lines-of– sight, the natural trend of the tissue structure has been preserved. To better visualize the effects of alignment, a small region containing blood flow from Figure 1 was aligned both within and between lines of sight and enlarged. Figure 2 shows the unaligned image, which contains a 40 mm vessel, and Figure 13 exhibits the aligned image. After alignment, blood flow is more easily distinguished from the surrounding tissue. Resulting velocity maps As previously mentioned, the mean velocity within each small region of a 2-D slice is estimated after alignment. Threshold values of the likelihood magnitude and postwall-filtered power are then applied, and the velocity is set to zero in regions in which either of these values is too low. Threshold levels are determined empirically to suppress noise-induced velocity values that fluctuate wildly over small depth ranges. A nonlinear filter is then employed to remove small regions of noise while pro-

Fig. 12. Aligned raw data of the ciliary body of a rabbit eye.

ducing smooth contours. Following this signal processing, the resulting velocity contours are overlaid onto the original aligned data set (Figure 12), and the regions containing blood flow are visualized, as shown in Figure 14. These overlaid contours are then evaluated to ensure that the identified regions of flow correspond to slanted RF data, which are indicative of blood flow, and to roughly estimate the blood velocity within the region based on the slant of the RF carrier. The peak axial velocity in Figure 14 is approximately 1 cm/s in the vessel on the left and 4 mm/s in the vessel on the right. The combination of alignment, high-resolution velocity estimation, and nonlinear postprocessing has allowed us to detect, resolve and maintain connectivity in these vessels within and between lines-of–sight. The image shown in Figure 14 is that produced by the averaging algorithm. Using the reference line algorithm, the results are similar; however, several extraneous regions of tissue motion trigger a low-velocity contour. Following the application of the incremental algorithm, significant regions of tissue motion with velocities of 1 mm/s or greater trigger false velocity contours. Thus, although each of these algorithms successfully maps axial velocities above several mm/s, maintaining continuity of small vessels requires the use of the refer-

822

Ultrasound in Medicine and Biology

Volume 24, Number 6, 1998

Fig. 13. Enlarged section of aligned ciliary body from Figure 12.

ence line or averaging alignment stategy, together with an algorithm that aligns structures between lines-of– sight. CONCLUSIONS Employing an ultrasound transducer center frequency of approximately 38 MHz, an effective spatial resolution of '40 mm and a velocity resolution of 1 mm/s can be reached. These very high spatial and velocity resolutions increase sensitivity to tissue motion, which can only be overcome by utilizing a preprocessing algorithm to align consecutive pulses with respect to each other. Three alignment algorithms were developed

and compared. To provide an absolute basis for performance comparison, we derived a confidence interval for the estimation of the correlation peak. Evaluation of this expression for our experimental conditions predicts a residual jitter after alignment on the order of 1 ns. Using specially constructed phantom data, our two measures of residual jitter (correlation magnitude and postwall-filtered power) demonstrate that correlation against a single reference line is superior to a nearestneighbor algorithm. By then expanding this algorithm to allow for averaging over multiple reference lines, alignment is further improved. In addition, this analysis demonstrates that data interpolation producing 1 ns pixels,

Ultrasonic mapping of microvasculature ● B. G. ZAGAR et al.

823

5 E$@p0~t 2 tr! 1 nr~t!# z @p0~t 2 ti 1 t! 1 ni~t!#% 5E

~t 2 t ! Ç ~t 2 t 1 t ! p Hp Ç J a b 0

r

i

0

5 E$ p 0~ a ! p 0~ a 1 t r 2 t i 1 t !% 5 R p0,p0~ t 2 @ t i 2 t r#!

(14)

with: E { . . . } expected value of { . . . }. Derivation of Equation (5): First, we develop an expression for the variance in the estimation of the peak magnitude of a time-limited correlation function. Second, the variance in the peak magnitude is used to determine the variance in estimation of the precise location of the correlation peak. Clearly, if the correlation function is very narrow (indicating the wide bandwidth of the process), an equal error in the magnitude will translate into a smaller error in the location as compared to a narrow band process. We assume the pulse, p 0 (t), has a power spectral density with bandwidth, B, centered at frequency, f 0 . The theoretical autocorrelation function for that power spectral density is given by: Fig. 14. Velocity profile contours overlaid onto raw aligned data in Figure 12.

which corresponds to our predicted residual jitter, allows for nearly optimal alignment, in that the postwall-filtered power cannot be further reduced. By applying these algorithms to an interpolated data set, the residual jitter in the image can be maintained well below 6 1/2 pixel from pulse to pulse within a single line-of–sight. Expanding the nearest-neighbor algorithm to allow for alignment between lines-of–sight further reduces tissue motion and allows for improved visualization of vessel continuity between lines-of–sight. Finally, we demonstrate that low-velocity blood flow can be reliably detected and mapped within closely spaced vessels. Even for flow velocities that are of equal magnitude to the velocity of physiological motion, the velocity estimator detected flow within the region of the vessel and did not detect motion in the surrounding tissue. APPENDIX To derive Equation (3): Rpr,pi~t! 5 E$pr~t!pi~t 1 t!%

(13)

R p0,p0~ t ! 5 AB

sin~2 p B t ! cos~2 p f 0t !, 2pBt

(15)

where A is a scaling constant to accommodate for the correct total power. Expanding this function into a Taylor-series around t 5 0 yields:

F

Rˆ p0,p0~ t ! 5 R p0,p0~0! 1 2

G

~ p B! 2 1 3~2 p f 0! 2 2 t . 6

(16)

Clearly, Rˆ p 0,p 0(0) is an unbiased estimate of R p 0,p 0(0) with a variance given by: var$Rˆp0,p0~0!% 5 E$@Rˆp0,p0~0! 2 Rp0,p0~0!#2%

HF

5 Rp20,p0~0! z E 1 2

F

5 Rp20,p0~0!

GJ

~pB!2 1 3~2pf0!2 2 t 21 6

G

~pB!2 1 3~2pf0!2 2 z E$t4%. 6

2

(17)

Assuming that the values of t follow a Gaussian distribution with zero mean and variance, s2(t), studies have shown (Bendat and Piersol 1986) that E{ t 4 } 5 3E { t 2 } 2 5 3 s 4 ( t ). Therefore, Equation (17) becomes:

824

Ultrasound in Medicine and Biology

var$Rˆ p0,p0~0!% 5 R p20,p0~0!

F

G

~ p B! 2 1 3~2 p f 0! 2 6

2

z 3 s 4~ t !. (18)

Rearranging terms yields:

s 2~ t ! 5 5

var$Rˆ p0,p0~0!% 6/ Î3 2z ~ p B! 1 3~2 p f 0! R p20,p0~0! 2

(19)

6/ Î3 z «$Rˆ p0,p0~0!%, ~ p B! 2 1 3~2 p f 0! 2

which is the expression in Equation (5). Acknowledgements—This work was supported in part by NIH Grant EY 11468, the Whitaker Foundation and the Austrian Ministry for Science and the Arts.

REFERENCES Baish JW, Gazit Y, Berk DA, Nozue M, Baxter LT, Jain RK. Role of tumor vascular architecture in nutrient and drug delivery: An invasion percolation-based network model. Microvasc Res 1996;51: 327–346. Bendat JS, Piersol AG. Random data: Analysis and measurement procedures. 2nd ed. New York John Wiley & Sons, Inc., 1986. Bigler SA, Deering RE, Brawer MK. Comparison of microscopic vascularity in benign and malignant prostate tissue. Human Pathol 1993;24:220 –226. Coleman DJ, Silverman RH, Rondeau MJ, Lizzi FL, McLean IW, Jakobiec FA. Correlations of acoustic tissue typing of malignant melanoma and histopathologic features as a predictor of death. Am J Ophthalmol 1990;110:380 –388. Coleman DJ, Silverman RH, Rondeau MJ, Coleman JA, Rosberger D, Ellsworth RM, Lizzi FL. Ultrasonic tissue characterization of uveal melanoma and prediction of patient survival after enucleation and brachytherapy. Am J Ophthalmol 1991;112:682– 688. Crochiere RE, Rabiner LR. Multirate digital signal processing. Englewood Cliffs, NJ: Prentice Hall Inc., 1983. Ferrara KW, Algazi VR. A new wideband spread target maximum likelihood estimator for blood velocity estimation, Part one: Theory. IEEE Trans Ultrason Ferro Freq Cont 1991a;38:1–16. Ferrara KW, Algazi VR. A new wideband spread target maximum likelihood estimator for blood velocity estimation, Part two: Evaluation of estimators with experimental data. IEEE Trans Ultrason Ferroelec Freq Cont 1991b;38:17–26.

Volume 24, Number 6, 1998 Ferrara KW, Zagar BG, Sokil-Melgar J, Algazi VR. High resolution 3D color flow mapping: Applied to the assessment of breast tumors. Ultrasound Med Biol 1996a;22:293–304. Ferrara KW, Zagar B, Sokil-Melgar J, Silverman R, Aslanidis Y. Estimation of blood velocity with high frequency ultrasound. IEEE Trans Ultrason Ferro Freq Cont 1996b;43:149 –157. Folberg R, Pe’er J, Gruman LM, Woolson RF, Jeng G, Montague PR, Moninger TO, Yi H, Moore KC. The morphologic characteristics of tumor blood vessels as a marker of tumor progression in primary human uveal melanoma: a matched case-control study. Human Pathol 1992;23:1298 –1305. Folberg R, Rummelt V, Parys-Van Ginderdeuren R, Hwang T, Woolson RF, Pe’e J, Gruman LM. The prognostic value of tumor blood vessel morphology in primary uveal melanoma. Ophthalmology 1993;100:1389 –98. Folkman J, Nerler E, Abernathy C, Williams G. Isolation of a tumor factor responsible for angiogenesis. J Exp Med 1971;33:275–280. Foltz RM, McLendon RE, Friedman HS, Dodge RK, Bigner DD, Dewhirst MW. A pial window model for the intracranial study of human glioma microvascular function. Neurosurgery 1995;36:976 – 985. Jensen JA. Estimation of blood velocities using ultrasound. Cambridge: Cambridge University Press, 1996. Kadi AP, Loupas T. On the performance of regression and stepinitialized IIR clutter filters for color Doppler systems in diagnostic medical ultrasound. IEEE Trans Ultrason Ferro Freq Cont 1995; 42:927–937. Less JR, Skalak TC, Sevick EM, Jain RK. Microvascular architecture in a mammary carcinoma: branching patterns and vessel dimensions. Cancer Res 1991;51:265–273. Less JR, Posner MC, Skalak TC, Wolmark N, Jain RK. Geometric resistance and microvascular network architecture of human colorectal carcinoma. Microcirculation 1997;4:25–33. Lockwood GR, Ryan LK, Hunt JW, Foster FS. Measurement of the ultrasonic properties of vascular tissues and blood from 35– 65 MHz. Ultrasound Med Biol 1991;17:653– 666. Marple SL. Digital spectral analysis with applications. Englewood Cliffs, NJ: Prentice Hall, Inc., 1987. Netti PA, Roberge S, Boucher Y, Baxter LT, Jain RK. Effect of transvascular fluid exchange on pressure-flow relationship in tumors: A proposed mechanism for tumor blood flow heterogeneity. Microvasc Res 1996;52:27– 46. Rummelt V, Folberg R, Rummelt C, Gruman LM, Hwang T, Woolson RF, Yi H, Naumann GO. Microcirculation architecture of melanocytic nevi and malignant melanoma of the ciliary body and choroid. A comparative histopathologic and ultrastructural study. Ophthalmology 1994;101:718 –27. Walker WF, Trahey GE. A fundamental limit on delay estimation using partially correlated speckle signals. IEEE Trans Ultrason Ferro Freq Cont 1995;42:301–308.