Ultrasound in Med. & BioL Vol. 18, No. 10, pp. 861-872, 1992 Printed in the U.S.A.
0301-5629/92 $5.00 + .00 © 1992 Pergamon Press Ltd.
OOriginal Contribution A METHOD FOR COMPUTER SIMULATION OF ULTRASOUND DOPPLER COLOR FLOW IMAGESnl. THEORY AND NUMERICAL METHOD A. T. KERR* a n d J. W. HUNT* tThe Ontario Cancer Institute and Department of Medical Biophysics, University of Toronto, 500 Sherbourne Street, Toronto, Ontario, Canada M4X 1K9, and *Ontario Cancer Treatment and Research Foundation, Kingston Regional Cancer Centre, King Street West, Kingston, Ontario, Canada, K7L 2V7
(Received 31 May 1991; in final form 6 July 1992) Abstract--Ultrasound imaging systems utilizing the pulsed Doppler principle are capable of providing images of blood flow in real time. We present a useful method for simulating flow images on a computer. Our method assumes that blood and surrounding tissue consist of many point-like scatterers positioned randomly in three dimensions. The position-dependent acoustic response of each scatterer is calculated using the acoustic impulse response method. This method takes into account the spatial effects of the transducer geometry on both the amplitude and temporal response of point-scattering. Details of theory, assumptions made in the simulation, and numerical methods are described fully for a spherically focused transducer, as well as a discussion of signal processing for generation of the flow image. Motion of a single scatterer is investigated to test the performance of the simulation algorithm. This simulation method could potentially be beneficial for detailed study of current and future flow imaging systems.
Key Words: Ultrasound, Doppler color flow imaging, Pulsed Doppler, Computer simulation.
INTRODUCTION
tion in velocity estimates (Angelson 1981). Despite this compromise, flow imaging provides valuable spatial information of velocity and direction of flow, and is gaining wide use in cardiology and peripheral vascular evaluation ( Kurjak et al. 1988; Merritt 1987; Roelandt 1986). The potential use of flow imaging in new areas, such as imaging of blood flow patterns in human tumors or monitoring transplant rejection, is being investigated (Adler et al. 1990; Cosgrove et al. 1990; Merritt 1987). To ensure that optimum systems are developed for these and possibly other applications, methods are needed that allow: 1) study of the current limitations and 2) testing of potential new improvements in flow imaging systems. A method widely used for studying B-scan imaging and standard pulsed Doppler ultrasound measurement is that of numerical simulation. Many investigators have presented simulated B-scan images (Bamber and Dickinson 1980; Foster et al. 1983), and from these images gained useful knowledge of the role of the transducer's geometry in the structural appearance of the image and in affecting image contrast. The use of simulation methods have been applied to study continuous-wave (CW) Doppler (Mo and Cobbold 1986a, 1986b), and pulsed Doppler (Azimi and
Real-time pulse-echo imaging, commonly known as B-scan imaging, allows dynamic visualization of softtissue anatomy. The desire to include information about flow dynamics (direction and velocity of flow) has led to the development of Doppler Color Flow Imaging (DCFI) (Bommer and Miller 1982; Kasai et al. 1985; Namekawa et al. 1982), referred to in this paper as flow imaging. Flow imaging is an adaptation of conventional multichannel pulsed Doppler techniques (Brandestini 1978 ), whereby the process of acquiring flow information along a single line of depth is extended to many lines in real time. The multichannel pulsed Doppler instrument (Brandestini 1978) sends successive bursts clown a single line of depth into tissue in order to obtain the required phase information for estimating velocity at a series of points along the line of depth. The flow imaging system does the same, except the process is repeated for each line of the image. Consequently, the requirement for realtime display means fewer successive bursts are sent down each line, leading to a greater statistical variaRequests for reprints should be sent to J. W. Hunt, The Ontario Cancer Institute, 500 Sherbourne Street, Toronto, Ontario, Canada, M4X IK9. 861
862
Ultrasound in Medicineand Biology
Kak 1985; Brown et al. 1985; Jorgensen et al. 1971 ). More recently, the case of flow imaging (Bonnefous and Pesque 1986) has been explored. Bonnefous and Pesque ( 1986 ) presented a simulation method based on a point-scattering model for testing quantitatively improved velocity estimator techniques for flow imaging applications. The velocity estimators were compared successfully with simulation of velocity profiles of axisymmetric flow. However, the model did not take into account the scanning nature of the transducer, and flow images were not simulated. We believe that in order to fully study flow imaging systems, the simulation method should be able to consider all the steps in flow image generation, and be able to produce actual flow images. We present in this paper a method for simulating the entire flow image making process. This approach considers a volume of many point-scatterers, or targets, positioned randomly in three dimensions. Blood is represented as moving targets, and vessel wall and surrounding tissue is represented by more strongly scattering stationary targets. The backscatter signal from the cloud of scatterers is calculated by summing together the responses from each of the many scatterers. The response from each target is determined in a similar way as described by Bonnefous and Pesque (1986). However, the scattering response is calculated using the method of acoustic impulse response (Arditi et al. 1981; Foster et al. 1983), rather than expressing the received scatter signal as the transmitted signal weighted by a position-dependent scaling factor. The method of acoustic impulse response takes into account both the change in amplitude and shape of the scatter signal as targets move from one position to another. This type of approach will be particularly important when studying strongly focused transducer geometries. The scanning of the transducer in the image plane and data acquisition are modeled realistically by keeping track of the position of all the scatterers (moving or stationary) while the transducer sends successive bursts down each line of depth, scans to a new position, and repeats the process until the image plane is covered. Applying the acoustic impulse response approach to flow image simulation introduces new and unique numerical challenges not encountered in conventional B-scan simulation. As will be described, several critical compromises are required in order that flow images can be simulated accurately within realistic computation and memory limits. Following a description of the simulation model and its numerical implementation on a computer, we present results of the treatment of a single moving scatterer in order to test the performance of our simulation.
Volume18, Number 10, 1992 THEORY AND UNDERLYING ASSUMPTIONS
Scattering nature of blood The predominant sources of scattering in blood are the red blood cells (RBCs) ( Reid et al. 1969 ). The average diameter of an RBC is 7 #m, which is small compared to the wavelength of clinical ultrasound (150-400 um). This implies that RBCs may act in scattering as a collection of independent random point-like scatterers. Based on theoretical predictions of backscatter power and comparison with measured values from blood, the point-scattering model has been shown to be valid for hematocrit below 25% (Shung et al. 1976). However, as shown by Shung et al. (1976), RBCs in blood begin to interact at a hematocrit above 25%, and in normal blood (45% hematocrit) RBCs are "strongly interacting." This interaction leads to a decrease in the scattering cross-section of blood. Studies of CW Doppler signals from normal blood under laminar flow conditions (Mo and Cobbold 1986a, 1986b; Mo 1991) indicate that a pointscattering model can still be used to model the statistical properties of scattered signals. However, the random point-scattering model does not predict correctly the total scattering power from normal blood. Therefore, the relative scattering power from blood and tissue must be corrected within the model based on values measured experimentally. The point-scattering model is also not expected to hold for turbulent flow conditions where it is known that both the scattering power and statistical nature of the scattered signals differ significantly from the result of random scattering (Bascom et al. 1988). The remaining assumptions about the scattering medium necessary for the development of our simulation method are listed as follows: 1. Since blood is a weak scatterer, multiple scattering is ignored. 2. The sound-propagating medium is uniform, with a constant speed of sound. 3. Attenuation effects due to absorption and scattering are neglected. 4. The relative amplitude of scattering between blood and the surrounding tissue is a constant. 5. The concentration of scatterers is invariant spatially. In assuming a uniform medium, the model is limited to representing situations where refraction effects at the vessel wall are not present. In addition, neglecting attenuation effects restricts the use of the model for phantom studies or studying vessels close to the surface. However, it will be shown in the description of the simulation approach that attenuation could in principle be included in the model. The theory of the acoustic impulse response ap-
Ultrasound Doppler color flowimaging--I • A. T. KERRand J. W. HUNT proach for calculating the response of a transducer when a point scatterer is in the field is described in the next section.
If a scatterer is placed at the field point ?, and assuming reciprocity of scattering (Arditi et al. 1981 ), Hunt et al. (1983) have reviewed that the output voltage of the transducer V(t) is expressed by
Acoustic impulse response This section describes the mathematical derivation of the backscatter signal from a cloud of scatterers. The time-varying pressure distribution p (7, t) in the field of a transducer (Arditi et al. 1981; Stephanishen 1971 ) can be expressed as
OVo(t)
p(~, t) = - o - -
Ot
.h(?, t),
(1)
where 7 is the position of the field point, o is the mass density of the sound-propagating medium, and t is time. The term, Vo(t) is the instantaneous normal particle velocity induced by an input excitation voltage Vi, (t) across the face of the transducer. The symbol • indicates a convolution. The function h(-~, t) is generally called the acoustic impulse response and is defined by
h(-?, t) = ~
1 f s 6(t-r'/c) r'
dS.
(2)
As shown in Fig. 1, r' is the distance from an incremental surface element d S to the field point at ?.
Spherical Radiator
Point / Scat'ererX';,
P,a.e ,y
Fig. 1. Geometry used for representing the acoustic impulse response and calculated backscatter signal for a stationary target. Arrangement of field points in the focal plane is for numerical simulation and is described in the text.
863
V(t) = K - o
OV~,(t) ] O----~*g'(t)*g2(t)*S(t) • ht(~, t).h~(7, t),
(3)
where g~(t) and g2(t) are the transmit and receive electromechanical impulse responses of the transducer system, respectively. The constant K has been added to the expression by Hunt et al. (1983) to take into account the different scattering strengths between moving blood and surrounding tissue. The term S(t) is the scattering impulse response which in our model is assumed to be position independent. The terms hi(7, t) and h,(7, t) are the transmit and receive acoustic impulse responses, respectively. Rather than calculate the terms in square brackets of eqn (3), we have modeled the pulse-echo response as a Gaussian-shaped pulse, which we represent as Pc(t). Therefore, eqn (3) can be rewritten as
V(t) = rPe(t)*h,(?, t ) . h f L t).
(4)
It should be noted that this simplification is not stringent and could be eliminated since for many transducer geometries the impulse response is a deltafunction at the natural focal point, and hence, Pe(t) could be measured experimentally. The useful property of the impulse response is that it separates the pulse generation stages of the system from the acoustic response. Therefore, any pulse shape can be analyzed in this way. Equation (4) describes the output voltage, also known as the backscatter or radio frequency (RF) signal from a single scatterer that is stationary. When the target velocity ~ is not zero, the Doppler effect leads to a slight frequency shift in the output voltage relative to that of the stationary scatter signal. However, as shown by Bonnefous and Pesque (1986), under most conditions 151 ~ c, so to a very good approximation the true response of a moving scatterer can be represented by the response of a static scatterer that undergoes the same displacement between transmission intervals as the moving scatterer. Extending eqn (4) to consider many scatterers, while assuming the same transducer is used on transmit and receive, yields N
V(t) = Pc(t)* Z Kiht(Ti, t)*hr(-?i, t),
(5)
i=l
where 7~ is the position of the ith scatterer, N is the
864
Ultrasoundin MedicineandBiology Volume18,Number10,1992
total number of scatterers, and Ki is the scale factor for the ith scatterer. Equation (5) describes the backscatter signal for the transducer held in one position, and for one excitation burst. This equation forms the basis of our computation algorithm, which further takes into account successive excitation bursts and transducer motion, to be described later in the paper. The transducer is modeled in this paper as a spherically focused source. In the next section, the numerical implementation of the impulse response is described which will be shown to require several simplifications of the transducer response. F O R M A T I O N OF D O P P L E R F L O W I M A G E S
Scattering volume geometry The geometry of the three-dimensional scattering volume as modeled for simulation of flow images is shown in cross-section in Fig. 2. The figure shows a cylindrical vessel containing moving scatterers, surrounded by a cloud of more strongly scattering stationary targets. The image plane in this case cuts the vessel longitudinally such that the vessel axis lies in the image plane. The scattering volume extends significantly far beyond the image plane so that contributions by targets away from the transducer axis are included ( Foster et al. 1983 ). The transducer is situated above and its direction of beam propagation makes an angle 0 with respect to the vessel axis. Scan
Direction
. •
•
/
Diameter
I Transducer
°• •.
•
[
"
•
°
• " • " "-
•
"
.
•
:. :
•
•
"
". " Static .
" "
°
*
:i"
:
•:
• "° ~: . ' '•f .
'='
°
.
-
•
•
°
..
. ~ . : ! ~ ' " .".": ".0:~"i~:."
The speed of the computation is inversely proportional to the number density of scatterers. To minimize the computation time, the m i n i m u m number of scatterers/unit volume is desired. In the simulation of B-scan images ( Foster et al. 1983), for transducers off/number ( focal length / aperture diameter) as low asf/1, there was found to be no difference in measurable image parameters for scatterer concentrations above 1 5 / m m 3. Normal blood contains approximately 5 × 106/mm 3 RBCs (Atkinson and Woodcock 1982), but the effective scatterer number density is unknown because many are in contact with each other. Consequently, as a conservative minimum, we have chosen to consider a number density for scatterers of about 2 0 / m m 3. The number of scatterers within the full-width-at-half-maximum ( F W H M ) sample volume is greater than 20, and the number of scatterers is greater than 200 within the - 4 0 dB sample volume for the transducer geometries considered, so Rayleigh scattering is achieved.
Transducer and scanning motion Flow images were simulated by modeling the transducer in the imaging system as a spherically focusing single element source• The primary reason for choosing this geometry is that it is a source for which the impulse response is known analytically (Arditi et al. 1981)• The scanning of the spherically focused transducer was modeled to be a linear translation• As shown in Fig. 2, the transducer moves laterally in a step-wise manner, and such that at each transducer position the beam axis maintains a constant angle 0 with respect to the vessel axis. When the transducer is modeled to move to the next line, it does so instantaneously• Hence, true mechanical scanning is not represented. Instead, the image format from this scanning motion is representative of linear array images. The simulation method could be modified to deal with mechanical or sector scan motion but is not critical to the purpose of this paper•
VelocitYProfile
Numerical approach •
•
•
i.'i"-.
/i I i
Depth Direction
Y
Fig. 2 . A two-dimensional representation of the three-dimensional scattering model. The spherically focused transducer scans linearly across the image plane. A cylindrical vessel is shown in cross-section containing moving scatterers at a fixed angle 0 with respect to the depth direction. The background scatterers are more strongly scattering and stationary.
According to eqn ( 5 ), the backscatter signal may be computed by summing all pulse-echo impulse response functions, corresponding to each target, before convolution with the pulse-echo response Pe(t)" This method has been successful and efficient for simulating B-scan images. The numerical approach of Foster et al. (1983) was to generate a look-up table (LUT) that stores sampled values of ht(~, t)*hr(?, t) for a finite series of field points arranged in a grid pattern (along the transducer axis and off-axis). A two-dimensional grid was feasible because the scattering properties of a spherical transducer are symmetric
Ultrasound Doppler color flow i m a g i n g - - I • A. T. KERR and J. W. HUNT
about the axis. Preliminary studies found that with our current memory and computation-time restraints, this method was not adequate for simulating flow images. The reason hinges on the pulsed Doppler signal being very sensitive to errors in the change of phase (at constant depth) of each scatter signal as the scatterers move in the acoustic field. The phase errors were primarily introduced due to the grid spacing, and due to the inability to represent the impulse response accurately in time. Grid spacing introduces the problem of spatial aliasing, and the finite sampling interval introduces the problem of temporal aliasing. To minimize phase errors, we have had to alter the LUT and the order of operation in eqn ( 5 ). The LUT for the present simulation stores the responses for a series of field points arranged along a single line in the plane z = 0 (Fig. 1). The coordinates for these field points are (x = 0, y = gAr, z = 0), g = 1, . . . . gmax, where Ar = 8 um is the spacing between field points, and gmx = 500 is the total number of field points. Thus, impulse responses were generated for points up to a maximum of 4 mm offthe central axis of the transducer. The use of a single line of points limits only the range above and below the line within which the true response of the acoustic source is represented. For each point on the line, a finely sampled discrete impulse response function hg(gAr, n a t ) is generated as follows: L-I
hg(gAr, n a t ) = ~, h(gAr, n a t + l A t / L ) , /=1 /7 = 1 . . . . .
Ng, g=
(6)
1,...,500,
where /xt is the sampling interval of the impulse response, and L is the number of samples within each /xt that are added to estimate the area under the analytical expression h(gAr, t) in the sampling interval At. The estimation of area between each sampling interval is required because convolution involves integration. The integer quantity Ng is the maximum number of samples required to represent the impulse response for the field point at y = gAr. The sampled impulse response is then convolved with itself and the pulse-echo response to yield Vg(gAr, n a t ) = hg(gAr, n a t ) • hg(gAr, n A t ) . P e ( n A t ) ,
n = 1. . . . .
Ng,
(7)
where Vg(gAr, n a t ) is the discrete representation of the resultant backscatter signal. The pulse-echo response band-pass filters the acoustic impulse response and allows Ve(gAr , n a t ) to be represented with a
865
larger sampling interval At', hence fewer number of values. The resampled scatter response is V'g(gAr, nat') = Vg(gAr, kAt), n = 1. . . . . Ng(At/At'), k = 1, /xt'//xt; . . . . Ng.
(8)
The values of L, At, and At' for all simulations were 500, 0.001 #s, and 0.05 #s, respectively. The Fourier transform of each discrete scatter signal is generated using a Fast-Fourier Transform (FFT) (Brigham 1974), and these values are stored in the LUT. Storing the frequency response values of the scatter signals allows the signals to be shifted exactly in time as the scatterer moves. For example, a scatterer at position ~ with coordinates (x, y, z) is represented by the LUT with a field point number g ( x , y) = NINT[ f ~ + yZ/Ar], where NINT is a nearest integer operator. The discrete frequency response of the scatter signal is then represented as A[-?, rnAo~] = A[g(x, y), mAw]e j~2z/c, m = 1. . . . .
1024, (9)
whereA[g(x, y), maw], m = 1. . . . . 1024, isthe FFT of V'g[g ( x , y ), nat'], n = 1. . . . . Ng( At / At'), and Aw is the frequency interval. The modulating term e j'~zz/c in the frequency domain corresponds to a time shift in the time domain. The time shift that is applied for the scatterer at ~ is 2 z / c where z is the z-coordinate of the scatterer. The main advantage to this type of numerical approach is that the grid spacing is very small and scatter signals are more correctly time shifted. Hence, phase errors in the pulsed Doppler signal are greatly reduced. The convenience of being able to study different pulse-echo responses without having to repeat the summation in eqn (5) is, however, lost. A flow chart (Fig. 3) shows how the algorithm proceeds once a LUT of frequency domain scatter response functions is generated. The algorithm was written in FORTRAN and ran on a fast computer (CRAY X-MP/24, Cray Research Corp, Eagan, MN, USA). A typical flow image requires about 10 min of central processing units (CPU) on the CRAY to complete. With use of a random number generator, the initial position of the ith target is assigned, where i = 1 . . . . . Nand Nis the total number of targets. A velocity of zero is assigned to targets if they lie outside the vessel boundary. The next step is to determine the position of the target (moving or stationary) relative to the transducer as the transducer's lateral position is stepped, and as it is modeled to send M successive bursts down each image line. For each successive
Ultrasound in Medicineand Biology
866
'Input Parameters: :~xb,+y~,+z ~ = target boundary NL = no. of image lines M = no. of bursts/line N --- no. of scatterers r, --- vessel radius
1 v,~ = maximum velocity n = flow profile constant I T = time between bursts. I 0 = Doppler angle W = width of image plane
'Random initial position of i'~ target: (x,,y~,z~)" Distance from target to vessel axis: r, = " 4 [ z ~ + ( y i / t a n O ) s i n O ] 2 + x i
' i=I,...N
2
If r,
0301-5629/92 $5.00 + .00 © 1992 Pergamon Press Ltd.
OOriginal Contribution A METHOD FOR COMPUTER SIMULATION OF ULTRASOUND DOPPLER COLOR FLOW IMAGESnl. THEORY AND NUMERICAL METHOD A. T. KERR* a n d J. W. HUNT* tThe Ontario Cancer Institute and Department of Medical Biophysics, University of Toronto, 500 Sherbourne Street, Toronto, Ontario, Canada M4X 1K9, and *Ontario Cancer Treatment and Research Foundation, Kingston Regional Cancer Centre, King Street West, Kingston, Ontario, Canada, K7L 2V7
(Received 31 May 1991; in final form 6 July 1992) Abstract--Ultrasound imaging systems utilizing the pulsed Doppler principle are capable of providing images of blood flow in real time. We present a useful method for simulating flow images on a computer. Our method assumes that blood and surrounding tissue consist of many point-like scatterers positioned randomly in three dimensions. The position-dependent acoustic response of each scatterer is calculated using the acoustic impulse response method. This method takes into account the spatial effects of the transducer geometry on both the amplitude and temporal response of point-scattering. Details of theory, assumptions made in the simulation, and numerical methods are described fully for a spherically focused transducer, as well as a discussion of signal processing for generation of the flow image. Motion of a single scatterer is investigated to test the performance of the simulation algorithm. This simulation method could potentially be beneficial for detailed study of current and future flow imaging systems.
Key Words: Ultrasound, Doppler color flow imaging, Pulsed Doppler, Computer simulation.
INTRODUCTION
tion in velocity estimates (Angelson 1981). Despite this compromise, flow imaging provides valuable spatial information of velocity and direction of flow, and is gaining wide use in cardiology and peripheral vascular evaluation ( Kurjak et al. 1988; Merritt 1987; Roelandt 1986). The potential use of flow imaging in new areas, such as imaging of blood flow patterns in human tumors or monitoring transplant rejection, is being investigated (Adler et al. 1990; Cosgrove et al. 1990; Merritt 1987). To ensure that optimum systems are developed for these and possibly other applications, methods are needed that allow: 1) study of the current limitations and 2) testing of potential new improvements in flow imaging systems. A method widely used for studying B-scan imaging and standard pulsed Doppler ultrasound measurement is that of numerical simulation. Many investigators have presented simulated B-scan images (Bamber and Dickinson 1980; Foster et al. 1983), and from these images gained useful knowledge of the role of the transducer's geometry in the structural appearance of the image and in affecting image contrast. The use of simulation methods have been applied to study continuous-wave (CW) Doppler (Mo and Cobbold 1986a, 1986b), and pulsed Doppler (Azimi and
Real-time pulse-echo imaging, commonly known as B-scan imaging, allows dynamic visualization of softtissue anatomy. The desire to include information about flow dynamics (direction and velocity of flow) has led to the development of Doppler Color Flow Imaging (DCFI) (Bommer and Miller 1982; Kasai et al. 1985; Namekawa et al. 1982), referred to in this paper as flow imaging. Flow imaging is an adaptation of conventional multichannel pulsed Doppler techniques (Brandestini 1978 ), whereby the process of acquiring flow information along a single line of depth is extended to many lines in real time. The multichannel pulsed Doppler instrument (Brandestini 1978) sends successive bursts clown a single line of depth into tissue in order to obtain the required phase information for estimating velocity at a series of points along the line of depth. The flow imaging system does the same, except the process is repeated for each line of the image. Consequently, the requirement for realtime display means fewer successive bursts are sent down each line, leading to a greater statistical variaRequests for reprints should be sent to J. W. Hunt, The Ontario Cancer Institute, 500 Sherbourne Street, Toronto, Ontario, Canada, M4X IK9. 861
862
Ultrasound in Medicineand Biology
Kak 1985; Brown et al. 1985; Jorgensen et al. 1971 ). More recently, the case of flow imaging (Bonnefous and Pesque 1986) has been explored. Bonnefous and Pesque ( 1986 ) presented a simulation method based on a point-scattering model for testing quantitatively improved velocity estimator techniques for flow imaging applications. The velocity estimators were compared successfully with simulation of velocity profiles of axisymmetric flow. However, the model did not take into account the scanning nature of the transducer, and flow images were not simulated. We believe that in order to fully study flow imaging systems, the simulation method should be able to consider all the steps in flow image generation, and be able to produce actual flow images. We present in this paper a method for simulating the entire flow image making process. This approach considers a volume of many point-scatterers, or targets, positioned randomly in three dimensions. Blood is represented as moving targets, and vessel wall and surrounding tissue is represented by more strongly scattering stationary targets. The backscatter signal from the cloud of scatterers is calculated by summing together the responses from each of the many scatterers. The response from each target is determined in a similar way as described by Bonnefous and Pesque (1986). However, the scattering response is calculated using the method of acoustic impulse response (Arditi et al. 1981; Foster et al. 1983), rather than expressing the received scatter signal as the transmitted signal weighted by a position-dependent scaling factor. The method of acoustic impulse response takes into account both the change in amplitude and shape of the scatter signal as targets move from one position to another. This type of approach will be particularly important when studying strongly focused transducer geometries. The scanning of the transducer in the image plane and data acquisition are modeled realistically by keeping track of the position of all the scatterers (moving or stationary) while the transducer sends successive bursts down each line of depth, scans to a new position, and repeats the process until the image plane is covered. Applying the acoustic impulse response approach to flow image simulation introduces new and unique numerical challenges not encountered in conventional B-scan simulation. As will be described, several critical compromises are required in order that flow images can be simulated accurately within realistic computation and memory limits. Following a description of the simulation model and its numerical implementation on a computer, we present results of the treatment of a single moving scatterer in order to test the performance of our simulation.
Volume18, Number 10, 1992 THEORY AND UNDERLYING ASSUMPTIONS
Scattering nature of blood The predominant sources of scattering in blood are the red blood cells (RBCs) ( Reid et al. 1969 ). The average diameter of an RBC is 7 #m, which is small compared to the wavelength of clinical ultrasound (150-400 um). This implies that RBCs may act in scattering as a collection of independent random point-like scatterers. Based on theoretical predictions of backscatter power and comparison with measured values from blood, the point-scattering model has been shown to be valid for hematocrit below 25% (Shung et al. 1976). However, as shown by Shung et al. (1976), RBCs in blood begin to interact at a hematocrit above 25%, and in normal blood (45% hematocrit) RBCs are "strongly interacting." This interaction leads to a decrease in the scattering cross-section of blood. Studies of CW Doppler signals from normal blood under laminar flow conditions (Mo and Cobbold 1986a, 1986b; Mo 1991) indicate that a pointscattering model can still be used to model the statistical properties of scattered signals. However, the random point-scattering model does not predict correctly the total scattering power from normal blood. Therefore, the relative scattering power from blood and tissue must be corrected within the model based on values measured experimentally. The point-scattering model is also not expected to hold for turbulent flow conditions where it is known that both the scattering power and statistical nature of the scattered signals differ significantly from the result of random scattering (Bascom et al. 1988). The remaining assumptions about the scattering medium necessary for the development of our simulation method are listed as follows: 1. Since blood is a weak scatterer, multiple scattering is ignored. 2. The sound-propagating medium is uniform, with a constant speed of sound. 3. Attenuation effects due to absorption and scattering are neglected. 4. The relative amplitude of scattering between blood and the surrounding tissue is a constant. 5. The concentration of scatterers is invariant spatially. In assuming a uniform medium, the model is limited to representing situations where refraction effects at the vessel wall are not present. In addition, neglecting attenuation effects restricts the use of the model for phantom studies or studying vessels close to the surface. However, it will be shown in the description of the simulation approach that attenuation could in principle be included in the model. The theory of the acoustic impulse response ap-
Ultrasound Doppler color flowimaging--I • A. T. KERRand J. W. HUNT proach for calculating the response of a transducer when a point scatterer is in the field is described in the next section.
If a scatterer is placed at the field point ?, and assuming reciprocity of scattering (Arditi et al. 1981 ), Hunt et al. (1983) have reviewed that the output voltage of the transducer V(t) is expressed by
Acoustic impulse response This section describes the mathematical derivation of the backscatter signal from a cloud of scatterers. The time-varying pressure distribution p (7, t) in the field of a transducer (Arditi et al. 1981; Stephanishen 1971 ) can be expressed as
OVo(t)
p(~, t) = - o - -
Ot
.h(?, t),
(1)
where 7 is the position of the field point, o is the mass density of the sound-propagating medium, and t is time. The term, Vo(t) is the instantaneous normal particle velocity induced by an input excitation voltage Vi, (t) across the face of the transducer. The symbol • indicates a convolution. The function h(-~, t) is generally called the acoustic impulse response and is defined by
h(-?, t) = ~
1 f s 6(t-r'/c) r'
dS.
(2)
As shown in Fig. 1, r' is the distance from an incremental surface element d S to the field point at ?.
Spherical Radiator
Point / Scat'ererX';,
P,a.e ,y
Fig. 1. Geometry used for representing the acoustic impulse response and calculated backscatter signal for a stationary target. Arrangement of field points in the focal plane is for numerical simulation and is described in the text.
863
V(t) = K - o
OV~,(t) ] O----~*g'(t)*g2(t)*S(t) • ht(~, t).h~(7, t),
(3)
where g~(t) and g2(t) are the transmit and receive electromechanical impulse responses of the transducer system, respectively. The constant K has been added to the expression by Hunt et al. (1983) to take into account the different scattering strengths between moving blood and surrounding tissue. The term S(t) is the scattering impulse response which in our model is assumed to be position independent. The terms hi(7, t) and h,(7, t) are the transmit and receive acoustic impulse responses, respectively. Rather than calculate the terms in square brackets of eqn (3), we have modeled the pulse-echo response as a Gaussian-shaped pulse, which we represent as Pc(t). Therefore, eqn (3) can be rewritten as
V(t) = rPe(t)*h,(?, t ) . h f L t).
(4)
It should be noted that this simplification is not stringent and could be eliminated since for many transducer geometries the impulse response is a deltafunction at the natural focal point, and hence, Pe(t) could be measured experimentally. The useful property of the impulse response is that it separates the pulse generation stages of the system from the acoustic response. Therefore, any pulse shape can be analyzed in this way. Equation (4) describes the output voltage, also known as the backscatter or radio frequency (RF) signal from a single scatterer that is stationary. When the target velocity ~ is not zero, the Doppler effect leads to a slight frequency shift in the output voltage relative to that of the stationary scatter signal. However, as shown by Bonnefous and Pesque (1986), under most conditions 151 ~ c, so to a very good approximation the true response of a moving scatterer can be represented by the response of a static scatterer that undergoes the same displacement between transmission intervals as the moving scatterer. Extending eqn (4) to consider many scatterers, while assuming the same transducer is used on transmit and receive, yields N
V(t) = Pc(t)* Z Kiht(Ti, t)*hr(-?i, t),
(5)
i=l
where 7~ is the position of the ith scatterer, N is the
864
Ultrasoundin MedicineandBiology Volume18,Number10,1992
total number of scatterers, and Ki is the scale factor for the ith scatterer. Equation (5) describes the backscatter signal for the transducer held in one position, and for one excitation burst. This equation forms the basis of our computation algorithm, which further takes into account successive excitation bursts and transducer motion, to be described later in the paper. The transducer is modeled in this paper as a spherically focused source. In the next section, the numerical implementation of the impulse response is described which will be shown to require several simplifications of the transducer response. F O R M A T I O N OF D O P P L E R F L O W I M A G E S
Scattering volume geometry The geometry of the three-dimensional scattering volume as modeled for simulation of flow images is shown in cross-section in Fig. 2. The figure shows a cylindrical vessel containing moving scatterers, surrounded by a cloud of more strongly scattering stationary targets. The image plane in this case cuts the vessel longitudinally such that the vessel axis lies in the image plane. The scattering volume extends significantly far beyond the image plane so that contributions by targets away from the transducer axis are included ( Foster et al. 1983 ). The transducer is situated above and its direction of beam propagation makes an angle 0 with respect to the vessel axis. Scan
Direction
. •
•
/
Diameter
I Transducer
°• •.
•
[
"
•
°
• " • " "-
•
"
.
•
:. :
•
•
"
". " Static .
" "
°
*
:i"
:
•:
• "° ~: . ' '•f .
'='
°
.
-
•
•
°
..
. ~ . : ! ~ ' " .".": ".0:~"i~:."
The speed of the computation is inversely proportional to the number density of scatterers. To minimize the computation time, the m i n i m u m number of scatterers/unit volume is desired. In the simulation of B-scan images ( Foster et al. 1983), for transducers off/number ( focal length / aperture diameter) as low asf/1, there was found to be no difference in measurable image parameters for scatterer concentrations above 1 5 / m m 3. Normal blood contains approximately 5 × 106/mm 3 RBCs (Atkinson and Woodcock 1982), but the effective scatterer number density is unknown because many are in contact with each other. Consequently, as a conservative minimum, we have chosen to consider a number density for scatterers of about 2 0 / m m 3. The number of scatterers within the full-width-at-half-maximum ( F W H M ) sample volume is greater than 20, and the number of scatterers is greater than 200 within the - 4 0 dB sample volume for the transducer geometries considered, so Rayleigh scattering is achieved.
Transducer and scanning motion Flow images were simulated by modeling the transducer in the imaging system as a spherically focusing single element source• The primary reason for choosing this geometry is that it is a source for which the impulse response is known analytically (Arditi et al. 1981)• The scanning of the spherically focused transducer was modeled to be a linear translation• As shown in Fig. 2, the transducer moves laterally in a step-wise manner, and such that at each transducer position the beam axis maintains a constant angle 0 with respect to the vessel axis. When the transducer is modeled to move to the next line, it does so instantaneously• Hence, true mechanical scanning is not represented. Instead, the image format from this scanning motion is representative of linear array images. The simulation method could be modified to deal with mechanical or sector scan motion but is not critical to the purpose of this paper•
VelocitYProfile
Numerical approach •
•
•
i.'i"-.
/i I i
Depth Direction
Y
Fig. 2 . A two-dimensional representation of the three-dimensional scattering model. The spherically focused transducer scans linearly across the image plane. A cylindrical vessel is shown in cross-section containing moving scatterers at a fixed angle 0 with respect to the depth direction. The background scatterers are more strongly scattering and stationary.
According to eqn ( 5 ), the backscatter signal may be computed by summing all pulse-echo impulse response functions, corresponding to each target, before convolution with the pulse-echo response Pe(t)" This method has been successful and efficient for simulating B-scan images. The numerical approach of Foster et al. (1983) was to generate a look-up table (LUT) that stores sampled values of ht(~, t)*hr(?, t) for a finite series of field points arranged in a grid pattern (along the transducer axis and off-axis). A two-dimensional grid was feasible because the scattering properties of a spherical transducer are symmetric
Ultrasound Doppler color flow i m a g i n g - - I • A. T. KERR and J. W. HUNT
about the axis. Preliminary studies found that with our current memory and computation-time restraints, this method was not adequate for simulating flow images. The reason hinges on the pulsed Doppler signal being very sensitive to errors in the change of phase (at constant depth) of each scatter signal as the scatterers move in the acoustic field. The phase errors were primarily introduced due to the grid spacing, and due to the inability to represent the impulse response accurately in time. Grid spacing introduces the problem of spatial aliasing, and the finite sampling interval introduces the problem of temporal aliasing. To minimize phase errors, we have had to alter the LUT and the order of operation in eqn ( 5 ). The LUT for the present simulation stores the responses for a series of field points arranged along a single line in the plane z = 0 (Fig. 1). The coordinates for these field points are (x = 0, y = gAr, z = 0), g = 1, . . . . gmax, where Ar = 8 um is the spacing between field points, and gmx = 500 is the total number of field points. Thus, impulse responses were generated for points up to a maximum of 4 mm offthe central axis of the transducer. The use of a single line of points limits only the range above and below the line within which the true response of the acoustic source is represented. For each point on the line, a finely sampled discrete impulse response function hg(gAr, n a t ) is generated as follows: L-I
hg(gAr, n a t ) = ~, h(gAr, n a t + l A t / L ) , /=1 /7 = 1 . . . . .
Ng, g=
(6)
1,...,500,
where /xt is the sampling interval of the impulse response, and L is the number of samples within each /xt that are added to estimate the area under the analytical expression h(gAr, t) in the sampling interval At. The estimation of area between each sampling interval is required because convolution involves integration. The integer quantity Ng is the maximum number of samples required to represent the impulse response for the field point at y = gAr. The sampled impulse response is then convolved with itself and the pulse-echo response to yield Vg(gAr, n a t ) = hg(gAr, n a t ) • hg(gAr, n A t ) . P e ( n A t ) ,
n = 1. . . . .
Ng,
(7)
where Vg(gAr, n a t ) is the discrete representation of the resultant backscatter signal. The pulse-echo response band-pass filters the acoustic impulse response and allows Ve(gAr , n a t ) to be represented with a
865
larger sampling interval At', hence fewer number of values. The resampled scatter response is V'g(gAr, nat') = Vg(gAr, kAt), n = 1. . . . . Ng(At/At'), k = 1, /xt'//xt; . . . . Ng.
(8)
The values of L, At, and At' for all simulations were 500, 0.001 #s, and 0.05 #s, respectively. The Fourier transform of each discrete scatter signal is generated using a Fast-Fourier Transform (FFT) (Brigham 1974), and these values are stored in the LUT. Storing the frequency response values of the scatter signals allows the signals to be shifted exactly in time as the scatterer moves. For example, a scatterer at position ~ with coordinates (x, y, z) is represented by the LUT with a field point number g ( x , y) = NINT[ f ~ + yZ/Ar], where NINT is a nearest integer operator. The discrete frequency response of the scatter signal is then represented as A[-?, rnAo~] = A[g(x, y), mAw]e j~2z/c, m = 1. . . . .
1024, (9)
whereA[g(x, y), maw], m = 1. . . . . 1024, isthe FFT of V'g[g ( x , y ), nat'], n = 1. . . . . Ng( At / At'), and Aw is the frequency interval. The modulating term e j'~zz/c in the frequency domain corresponds to a time shift in the time domain. The time shift that is applied for the scatterer at ~ is 2 z / c where z is the z-coordinate of the scatterer. The main advantage to this type of numerical approach is that the grid spacing is very small and scatter signals are more correctly time shifted. Hence, phase errors in the pulsed Doppler signal are greatly reduced. The convenience of being able to study different pulse-echo responses without having to repeat the summation in eqn (5) is, however, lost. A flow chart (Fig. 3) shows how the algorithm proceeds once a LUT of frequency domain scatter response functions is generated. The algorithm was written in FORTRAN and ran on a fast computer (CRAY X-MP/24, Cray Research Corp, Eagan, MN, USA). A typical flow image requires about 10 min of central processing units (CPU) on the CRAY to complete. With use of a random number generator, the initial position of the ith target is assigned, where i = 1 . . . . . Nand Nis the total number of targets. A velocity of zero is assigned to targets if they lie outside the vessel boundary. The next step is to determine the position of the target (moving or stationary) relative to the transducer as the transducer's lateral position is stepped, and as it is modeled to send M successive bursts down each image line. For each successive
Ultrasound in Medicineand Biology
866
'Input Parameters: :~xb,+y~,+z ~ = target boundary NL = no. of image lines M = no. of bursts/line N --- no. of scatterers r, --- vessel radius
1 v,~ = maximum velocity n = flow profile constant I T = time between bursts. I 0 = Doppler angle W = width of image plane
'Random initial position of i'~ target: (x,,y~,z~)" Distance from target to vessel axis: r, = " 4 [ z ~ + ( y i / t a n O ) s i n O ] 2 + x i
' i=I,...N
2
If r,
