FAST AND ACCURATE TRANSIENT ULTRASOUND PROPAGATION AND B-MODE IMAGING SIMULATION METHODS By Yi Zhu A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical Engineering – Master of Science 2013 ABSTRACT FAST AND ACCURATE TRANSIENT ULTRASOUND PROPAGATION AND B-MODE IMAGING SIMULATION METHODS By Yi Zhu Computer simulation plays a significant role in the development of diagnostic ultrasound sound systems. Traditional calculation methods for continuous wave ultrasound calculations converge slowly when applied to transient wave propagation. Fast and accurate calculation methods are developed and discussed in this thesis to improve the performance of computing transient ultrasound wave propagation. B-mode imaging is widely used in practical diagnostic ultrasound applications. The imaging system is based on the pulse echo model. Transient ultrasound pulse waves are generated from the transducer, and echo signals are received to generate images. Traditional imaging simulation methods use the impulse response method, which requires high sampling frequencies to reach acceptable accuracies. By employing a convolutional simulation model and using fast and accurate transient ultrasound wave calculation methods, the simulation performance can be improved significantly both in image quality and simulation speed. Copyright by YI ZHU 2013 ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Robert McGough, for his excellent guidance and support during my master’s degree, and for providing a nice atmosphere to study and conduct research in the Biomedical Ultrasonics and Electromagnetics Lab. Without his knowledge and insights in biomedical ultrasound, signal processing, ultrasound propagation and imaging theories, this thesis would not have been possible. I would like to thank Dr. Lalita Udpa for sharing her knowledge and experience in the imaging processing, as well as for proposing inspiring questions to my research. I would like to thank Dr. Selin Aviyente, whose persistent help and concerns were invaluable to this thesis. I would also like to thank my lab mates for their help in this thesis. Thanks Peter Bread for providing great revising advises. Many thanks to Pedro Nariyoshi, Yiqun Yang and Xiaofeng Zhao, without whom it would be a lonely lab. Finally, I would like to thank my parents, who were always supportive to my decisions. Especially thanks to my mom, for being accompany with me during the hardest time. iv TABLE OF CONTENTS LIST OF TABLES ......................................................................................................................... vii LIST OF FIGURES ....................................................................................................................... viii Chapter 1 Introduction ..............................................................................................................1 1.1 Motivation ..............................................................................................................1 1.2 Background...........................................................................................................2 1.2.1 Time-harmonic wave calculations ......................................................................2 1.2.2 Transient wave calculations ................................................................................2 1.2.3 B-mode imaging simulations ...............................................................................4 1.2.4 Ultrasound simulators ..........................................................................................5 1.3 Thesis Organization .............................................................................................5 Chapter 2 Theory ......................................................................................................................7 2.1 Transient wave propagation models ..................................................................7 2.1.1 Impulse Response ...............................................................................................7 2.1.2 Fast Nearfield Method .......................................................................................13 2.1.3 Time-Space Decomposition ..............................................................................17 2.1.4 Frequency Domain Time-Space Decomposition ............................................20 2.2 B-mode imaging simulation ...............................................................................23 2.2.1 Impulse Response Pulse-echo Model..............................................................23 2.2.2 Algorithms for B-mode imaging simulations ....................................................24 2.3 Error Evaluation ..................................................................................................27 Chapter 3 Simulation Results ................................................................................................29 3.1 Transient wave propagation simulations .........................................................29 3.1.1 Simulation Parameters ......................................................................................29 3.1.2 Reference pressure field ...................................................................................31 3.1.3 Impulse Response .............................................................................................32 3.1.4 Time-Space Decomposition ..............................................................................34 3.1.5 Frequency Domain Time-Space Decomposition ............................................36 3.1.5.1 Spectral Clipping .................................................................................37 3.1.5.2 Calculation error – time ......................................................................39 3.1.6 Comparison ........................................................................................................41 3.2 B-mode imaging simulations .............................................................................43 3.2.1 Simulation Parameters ......................................................................................44 3.2.2 Point target scattering ........................................................................................46 3.2.2.1 Single focal zone for transmit and receive apertures ......................46 v 3.2.3 3.2.2.1.A.. B-mode imaging simulation in FOCUS without system clock 46 3.2.2.1.B.................... Implementing algorithms with 40MHz system clock 52 3.2.2.1.C .....................Influence of sampling frequency and system clock 55 3.2.2.2 Multiple focal zones for transmit and receive aperture ....................60 Cyst Phantom scattering ...................................................................................66 Chapter 4 Discussions ...........................................................................................................70 4.1 Transient wave propagation models ................................................................70 4.1.1 Impulse Response .............................................................................................70 4.1.2 Fast Nearfield Method .......................................................................................71 4.1.3 Time-Space Decomposition ..............................................................................71 4.1.4 Frequency Domain Time-Space Decomposition ............................................72 4.2 B-mode imaging simulations .............................................................................72 4.2.1 Point target scattering ........................................................................................73 4.2.1.1 Single Focal Zone ...............................................................................73 4.2.1.2 Multiple focal zones ............................................................................74 4.2.2 Cyst Phantom scattering ...................................................................................74 Chapter 5 Conclusions ...........................................................................................................76 Chapter 6 Future Work...........................................................................................................77 REFERENCES .......................................................................................................................78 vi LIST OF TABLES Table 2-1 Basis functions of TSD for Hanning weighted signal ..................................... 19 Table 3-1 Comparison of calculation time .................................................................... 43 Table 3-2 Sampling frequency and system clock effect ................................................ 60 Table 3-3 Point target scattering B-mode imaging simulation comparison .................... 65 vii LIST OF FIGURES Figure 2-1 Rectangular Piston and the projections of observation points (a) – (e) used in impulse response calculation. ...............................................................................................11 Figure 2-2 Impulse responses calculated for a rectangular piston with b/a = 1.5 at 5 different locations that share a common z coordinate with z/a = 10. ................................12 Figure 2-3 Aperture Movement. ............................................................................................25 Figure 3-1 Rectangular transducer piston and simulation coordinates. ............................30 Figure 3-2 Reference Pressure Field. ..................................................................................32 Figure 3-3 Impulse Response Calculation Time vs. Error. .................................................34 Figure 3-4 Time-Space Decomposition Calculation Time vs. Error. .................................36 Figure 3-5 Spectral Clipping. ................................................................................................38 Figure 3-6 Frequency Domain Time-Space Decomposition Calculation Time vs. Error. 40 Figure 3-7 Comparison of all propagation models. .............................................................42 Figure 3-8 Reference B-mode simulated image for single focal zone simulation without system clock. ..........................................................................................................................47 Figure 3-9 Field II and FOCUS simulated image comparison at 25MHz sampling frequency.................................................................................................................................48 Figure 3-10 Field II and FOCUS simulated image comparison at 50MHz sampling frequency.................................................................................................................................49 Figure 3-11 Field II and FOCUS simulated image comparison at 100MHz sampling frequency.................................................................................................................................51 Figure 3-12 Reference B-mode simulated image for single focal zone simulation with 40MHz system clock. .............................................................................................................53 Figure 3-13 Comparison of implementing algorithms. .......................................................54 Figure 3-14 Comparison of simulated images generated with different sampling frequencies..............................................................................................................................57 Figure 3-15 Comparison of simulated images generated with different system clocks. 59 viii Figure 3-16 Reference B-mode simulated image for 5 transmit and receive zones........61 Figure 3-17 Comparison of Field II and third algorithm at 40MHz sampling frequency with 40MHz system clock. .....................................................................................................63 Figure 3-18 Comparison of Field II and third algorithm at 200MHz sampling frequency with 40MHz sampling frequency. ..........................................................................................64 Figure 3-19 Cyst phantom simulated image generated by Field II. ..................................67 Figure 3-20 Cyst phantom simulated image generated by the third algorithm (ALG3). ..68 Figure 3-21 Cyst Phantom RF data comparison. ................................................................69 ix Chapter 1 1.1 Introduction Motivation Computer simulations are often used for evaluating the performance of ultrasound models, beamformers, and signal processing methods as applied to medical ultrasound systems. As new ultrasound systems are created, the amount of RF data to be processed in computer simulations is increasing. The simulation time and memory usage expand accordingly as the size of temporal and spatial grids increase. For larger and more complicated problems, fast, accurate and memory-efficient ultrasound simulation models and algorithms are needed. In therapeutic ultrasound, calculations of time-harmonic ultrasound waves are important for computing the power deposition and subsequent heat transfer calculations. Traditional calculation models such as the Rayleigh-Sommerfeld integral are widely used in high intensity focused ultrasound (HIFU) applications [1]. When applied to transient pressure field calculations, which are important for ultrasound imaging and some therapeutic ultrasound applications, the slow convergence in the nearfield region makes the Rayleigh-Sommerfeld integral a less ideal calculation model [2]. Fast and accurate calculation models for transient pressure field are desired for therapeutic ultrasound applications. Efficient simulation models and algorithms are also desired for B-mode ultrasound imaging employing fast transient pressure calculation methods. 1 1.2 Background 1.2.1 Time-harmonic wave calculations Time-harmonic ultrasound waves are calculated for therapeutic ultrasound applications. In homogeneous media, the conventional calculation model is the Rayleigh-Sommerfeld integral [3][4]. Due to the slow convergence rate of the Rayleigh-Sommerfeld integral, equivalent integral approaches such as the rectangular radiator method [5] and the spatial impulse response method [6] have been developed. Closed-form spatial impulse response functions have been derived for common transducer geometries such as circular pistons [7], rectangular pistons [6], triangular pistons [8], and spherical shells [9][10][11]. For 2D and 3D problems and larger computational grids, the angular spectrum approach has been developed for high speed calculations [12][13]. The Fast Nearfield Method (FNM) speeds up the calculations in the nearfield region by combing integrands and reducing the number of integrals [2][15]. In an inhomogeneous medium, a hybrid angular spectrum method has been developed for large grid calculations [14]. Finite difference methods [16] have also been developed for nonlinear time-harmonic wave calculations. 1.2.2 Transient wave calculations In ultrasound imaging, short pulses are generated by an ultrasound transducer [17]. Transient pressure calculations are therefore important for ultrasound imaging simulations. The impulse response method is a popular approach for transient pressure field calculations [7][18][19][20]. In order to obtain the transient pressure field, a single convolution is executed between the time derivative of the input velocity signal and the spatial impulse response function for a specific transducer geometry. To extend the 2 impulse response to arbitrary transducer geometries, an approach of dividing the transducer surface into sub-rectangle or sub-triangles and summing up the response has been developed [21][22]. Due to the discontinuities at the edge of the transducer, the impulse response method often requires high a sampling frequency to avoid aliasing problems and to achieve an acceptable accuracy. The performance of the impulse response method is improved with advanced approaches [23][24]. Even better performance is achieved when Time-Space Decomposition (TSD) [25][26] is combined with the Fast Nearfield Method for transient pressure field calculations. Since the Fast Nearfield Method removes the singularities from each the integrals, the Fast Nearfield Method naturally avoids the aliasing problems that occur in the impulse response method and therefore achieves high accuracies with low sampling frequency requirements. Time-Space Decomposition, which divides the velocity signal into separate temporal and spatial functions, further improves the performance of the Fast Nearfield Method when calculating transient pressure fields. In standard Time-Space Decomposition, explicit temporal and spatial functions are available only for limited signal types. Frequency Domain Time-Space Decomposition (FDTSD) extends this approach to enable transient calculations with the Fast Nearfield Method for arbitrary input signal waveforms [27]. After the Discrete Fourier Transform (DFT)is applied to the input signal waveform, the signal is expressed by several discrete frequency domain components. Spectral clipping can control the trade-off between calculation speed and accuracy in Frequency Domain Time-Space Decomposition. Previously, the FDTSD method was developed for circular transducers. In this thesis, the FDTSD will be extended to rectangular transducers. 3 1.2.3 B-mode imaging simulations Due to the relatively low cost, the ability to perform real-time imaging, and the absence of ionizing radiation, ultrasound imaging has certain advantages over other imaging technologies [28]. Diagnostic ultrasound imaging modes include B-mode (brightness mode), C-mode, M-mode (motion mode) and several Doppler modes [17][29]. B-mode imaging is of particular interest for this thesis. In B-mode imaging, short pulses generated by piezo-electric transducers propagate in soft tissue. Echo signals are received by the same transducers and then processed to compose B-mode images. An impulse response model has been proposed for B-mode imaging [29]. The total impulse response of the system can be calculated by convolving individual impulse responses generated by different transducers, accounting for diffraction and scattering. For point targets, Jenson introduced a simplified impulse response model [30]. An equivalent pulse-echo signal can be obtained by convolving the second derivative of the excitation signal with the spatial transmit and receive impulse responses [31]. If the excitation signal can be expressed as a convolution between two identical signals, the simplified impulse response model can be then expressed in terms of transmit and receive pressure signals. The Fast Nearfield Method combined with Time-Spatial Decomposition/Frequency Domain Time-Spatial Decomposition is then applied to the simulation model to accelerate the simulation and reduce the nearfield error. Several Bmode imaging simulation algorithms will be developed in this thesis, and these will be compared to an impulse response model and to Field II [32]. 4 1.2.4 Ultrasound simulators Field II (http://field-ii.dk/) is an ultrasound simulation program [32] based on analytical expressions for the impulse response derived by Tupholme[33] and Stepanishen[18]. In the Field II program, the calculation is accelerated by dividing the transducer surface into rectangles and applying the far-field approximation. FOCUS (http://www.egr.msu.edu/~fultras-web/) is a cross-platform ultrasound simulation tool for rapidly and accurately calculating pressure fields. In FOCUS, calculations for continuous wave pressure in a large 3-D grid are further accelerated by the Angular Spectrum Approach, and for transient pressures the simulations are accelerated by the Fast Nearfield Method combined with Time-Space Decomposition [25] or combined with Frequency Domain Time-Space Decomposition [27]. The need for fast and accurate B-mode image simulation routines in FOCUS motives the research presented in this thesis. 1.3 Thesis Organization This thesis consists of 6 chapters. Chapter 1 presents some background on diagnostic ultrasound and explains the motivation for the research in this thesis. Chapter 2 gives basic theories of propagation and B-mode image simulation models. Several important formulas, derivations, and concepts are explained carefully in this chapter. For transient pressure field calculations, the impulse response, Fast Nearfield Method, Time-Space Decomposition, and Frequency Domain Time-Space Decomposition will be presented. For B-mode image simulations, the pulse-echo system simulation model, the convolutional waveform expression, and several different analytically equivalent algorithms will be described. This chapter also gives the error evaluation formula used 5 for evaluations of the simulation results. In chapter 3, simulations are generated based on formulas and algorithms presented in chapter 2. Simulation results are illustrated with figures, tables, and careful explanations. Chapter 4 discusses the simulation results and evaluates comparisons between different models and algorithms. Conclusions are drawn in Chapter 5. In chapter 6, some plans for the future work are described. 6 Chapter 2 2.1 Theory Transient wave propagation models In this section, the impulse response, Fast Nearfield Method, Time-Space Decomposition and Frequency Domain Time-Space Decomposition will be carefully explained. The rectangular transducer is the transducer shape source described in each sub-section. 2.1.1 Impulse Response The impulse response method can be derived from the Green’s function [18]. For a time-dependent pressure generated in a half-plane due to a vibrating piston on an infinite planar baffle, the velocity potential ( where ) ( ∫ ∫ indicates the source points and area of piston surface, ( ( ) obtained from the Green’s function is ) ( (2.1) ) indicates the observation points, ) is the piston surface velocity, and ( is the ) is the Green’s function. For the region outside the piston surface on the infinite baffle, the normal velocity is zero. The medium in which the sound wave is propagating is assumed to be isotropic with a constant propagation velocity. The initial pressure and derivative of the pressure is equal to zero. The pressure is then obtained from the momentum equation, ( where ) ( ) is the density of the homogenous medium. 7 (2.2) For uniform piston surface velocity ( ), equation (2.1) is equivalent to ( ) ( ) ( ∫ ∫ (2.3) ) The Green’s function for this problem is ( where ( ) ) (2.4) is the speed of sound in the propagation medium, and distance from the source point to the observation point. is the in the denominator represents the radiation into the half-plan. The velocity potential can be then expressed as ( ) ( ( ) ∫ ∫ ) (2.5) If the spatial integration in equation (2.5) is evaluated first, the integral can be described in a convolutional form: ( ) ( ) ( ) (2.6) where ( ) ∫ ( ) (2.7) is defined as the impulse response function of the piston evaluated at a single observation point. To obtain the pressure field at the observation point, equations (2.2) and (2.6) are combined, which yields 8 ( ) ̇( ) ( ) (2.8) The pressure field in the half-plane can now be obtained by evaluating the convolution between the time derivative of the surface velocity and the impulse response function. The convolution can be calculated numerically with the Fast Fourier Transform (FFT). Lockwood and Willette [6] gives a closed-form solution for the impulse response of a rectangular piston. An observation point ( point ) in the pressure field has a projection ) on the piston surface plane. Assume the rectangular piston is 2b in ( length and 2a in height. The piston surface plane is divided into four rectangles, where each has a corner at . The pressure field is obtained by combining the contributions of four sub-rectangles. If is outside of the piston region, a larger rectangle is defined that includes the projection point and the contributions of other rectangles that are subtracted in the calculation. After defining the longer sides of th sub-rectangle as and shorter sides as , the exact solution for the impulse response of the rectangular piston can be written as ( ) ∑ ( ) {∫ ∫ { ∫ { √ ( ) √ ( ) where 9 (2.9) } ( ) } ( ) } √ (2.10) √ √ { Whether a term is added or subtracted in equation (2.9) is determined by the location of the observation point relative to the rectangular piston. If the projection point is located at the edge of the rectangular piston, the contributions of the two sub-rectangles with zero shorter sides are equal to zero. Figure 2-1 shows the projections of observation points (a) - (e) on the piston surface plane. For each observation point, the pressure field is the combinations of contributions of four sub-rectangles. In Figure 2-1, region I – IV are defined for point (c). The width of the rectangular piston is and the height is share a common z coordinate with . All of the observation points . Point (a) is located outside of the paraxial region at x = 1.4a, y = 1.25a. Point (b) is located at the edge of the paraxial region at x = a, y = 0.75b. Point (c) is located inside of the paraxial region at x = 0.6a, y = 0.75b. Point (d) is located on the x-axis at x = 0.6a, y=0. Point (e) is located at x = 0, y = 0. The corresponding impulse response functions are shown in Figure 2-2. 10 Figure 2-1 Rectangular Piston and the projections of observation points (a) – (e) used in impulse response calculation. In the figure, a is half the width of the transducer, and b is half the height of the transducer. Region I – IV are defined for point (c). All of the observation points share a common z coordinate with . For point (a), x = 1.4a, y = 1.25a. For point (b), x = a, y = 0.75b. For point (c), x = 0.6a, y = 0.75b. For point (d), x = 0.6a, y=0. For point (e), x = 0, y = 0. 11 Figure 2-2 Impulse responses calculated for a rectangular piston with b/a = 1.5 at 5 different locations that share a common z coordinate with z/a = 10. The horizontal axis represents the time in micro second, and the vertical axis represents the 12 Figure 2-2 (cont’d) amplitude of the impulse response. The impulse responses in subfigures (a) – (e) are corresponding to observation points (a) – (e) in Figure 2-1. The main disadvantage of impulse response calculations is that aliasing limits the accuracy when the sampling frequency is low. For example, discontinuities in the impulse response produce large errors at low sampling frequencies. To calculate the pressure with the impulse response method, the impulse response function needs to be evaluated during the entire duration of the impulse response signal. The computation time and memory usage can be very large when the sampling frequency is high, which is needed to obtain a more accurate result. In Field II, the pressure field is calculated based on the impulse response evaluated in the far-field region. Hence, the calculation accuracy is also limited by the sampling frequency in Field II. 2.1.2 Fast Nearfield Method Lockwood and Willette [6] defines the steady-state acoustic field produced by a transducer piston for a time-harmonic excitation as ( where ) ( ) (2.11) is the center frequency of the excitation in radians per second, density of medium, is the constant normal velocity on the piston surface, coordinate of the observation field point, and ( is the is the ) is the Fourier transform of the impulse response. For a rectangular piston [15], the Fourier transform of the impulse response over the corner of a rectangle can be written as 13 ( ) ( √ [ ] √ ∫ { √ } √ (2.12) √ ∫ { √ where variables } √ ) are the shorter and longer sides of the rectangle, respectively. After changing and replacing the wavenumber with , where is the speed of sound, the expression (2.12) is equivalent to the solution given by Lockwood and Willette [6]. The total Fourier transform of the impulse response for a rectangular piston is ( where ) ( ∑ ) (2.13) are the shorter and longer sides defined for up to four sub-rectangles. The sign of each contribution is decided by whether the corresponding sub-rectangle is added or subtracted. The Fast Nearfield Method for a rectangular transducer is obtained after three improvements are applied to equation (2.12). The first improvement involves deriving equivalent integrals from (2.12) and subtracting the singularities. The second improvement reduces repeated calculations in the new integral expressions. By combining integrals that share same integrand and common upper or lower limits, the 14 third improvement is realized. After combining these three steps, the pressure field is evaluated with Gauss quadrature. The expression that describes the Fast Nearfield Method for pressures evaluated over one corner of a rectangular piston is given by √ ( ) ( ∫ (2.14) √ ∫ ) Combining the shared integrand of sub-rectangles, the Fast Nearfield Method gives the expression for calculating pressure field generated by a rectangular piston as ( ) {( √ )∫ ( ) ( ) )∫ ( ) )∫ ( ) ( )∫ ( ) 15 ( ) ( ( ) √ ( ) ( ) √ ( ) ( √ ( ( ( ) ( ) ( ( (2.15) ) ) ) ) } In Equation (2.15), is half the width and is half the height of the rectangular piston. The sign of the contribution of each sub-rectangle is automatically determined by leading terms. The combination of terms with the same integrands reduces the number of integrals from eight to four in nearfield calculations with this expression. For points located on the edge of the transducer, the contributions from the two rectangles with zero-length short edges are equal to zero, and only two remaining integrals need to be evaluated. Equation (2.15) also subtracts singularities to improve the rate convergence, which is much better than the impulse response method. For a transient pulse excitation ( ), the Fast Nearfield Method is given as ( ) {∫ ∫ ( ) ( ) ( ( ) ) ( ( ) ) ( ( ) ) (2.16) ∫ ∫ ( ( ) ( ) ( ) ( ) where 16 ) ( ) ( ) ( ) } √( ) ( ) √( ) ( ) √( ) ( ) √( ) ( ) (2.17) { In Equation (2.16), the calculation involves both the time variable and variables that contain only spatial information. To further improve speed of the computation, the TimeSpace Decomposition method is combined with the Fast Nearfield Method. 2.1.3 Time-Space Decomposition Time-Space Decomposition converts the integrals in the Fast Nearfield Method, where both temporal and spatial variables are involved, into an analytical equivalent sum of spatial integrals that are weighted by time-dependent factors [25]. Since the time dependence is removed from the integrals, the equivalent superposition of spatial integrals converges rapidly. To further accelerate the calculation, the repeated terms in the integrals are calculated once, stored, and re-used in the following calculations. The resulting combination of the Fast Nearfield Method with Time-Space Decomposition is therefore fast and accurate in the nearfield region. For a temporally windeowed transient pulse, Time-Space Decomposition can be applied to the function ( ): ( ) ( )∑ 17 ( ) ( ) (2.18) where if ( ) , and otherwise. The variable ( ) only contains spatial information of the observation point coordinates, thus the function ( be evaluated by the temporal functions ( ) and spatial functions ) can ( ) seperately. By applying Equation (2.18), Equation (2.16) can be expressed as ( ) {∑ ( ) (∫ ( ∫ ( ) ( ) ) ( ) ) (2.19) ( ) (∫ ∫ )} where ( ) - ( ) ( ) and ) ( ) { ( ( ) ( ) ( ) are the same in Equation (2.17). The closed-form expressions of (2.20) ( ) and ( ) for Hanning-weighted pulse are given by D. Chen and R. J. McGough [26]. For 18 Hanning-weighted pulse, . The temporal and spatial functions for a Hanning- weighted pulse [26] are listed in Table 2-1. Table 2-1 Basis functions of TSD for Hanning weighted signal Temporal basis function ( ) ( ( ) ( ( ) ( ( ) ) ( ( ) ( ( ) ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ) ) ( ) ) ( ( ) ) ( ( ( ) ) ) ( ) ( ) ) ( ( ) Spatial functions ( ) ( ) ( ) When calculating nearfield pressures with the Fast Nearfield Method combined with Time-Space Decomposition, the temporal basis functions stored for all time values. The variable that contains the spatial information is then evaluated for each coordinate. The spatial functions after each ( ) are evaluated first and is calculated. The expression ( ( ) are evaluated immediately ) can also be evaluated directly with the current time point and observation coordinates by the direct evaluation of expression for the input signal. After all and ( ( ) are evaluated, the integrals containing ) are calculated by Gauss quadrature. Time-Space Decomposition then separates the transient pulse signal, combines integrals with the same integrands, 19 and exploits repeated calculations. This combination speeds up the Fast Nearfield Method without adversely impacting the accuracy of the nearfield pressure calculation. 2.1.4 Frequency Domain Time-Space Decomposition When combined with Time-Space Decomposition, the Fast Nearfield Method generates faster and more accurate transient pressure results in the nearfield region. The limitation of standard Time-Space Decomposition is that exact expressions for and ( ) ( ) are obtained only for a few signal types. For exponential functions containing non-integer powers of ( ) or if the power is other than one, it is impossible to expand the function ( ) into separate temporal and spatial functions. In order to exploit the advantage of Time-Space Decomposition for arbitrary single representations, Frequency Domain Time-Space Decomposition has been developed. Frequency Domain Time-Space Decomposition method is derived for circular piston by E. J. Alles et al [27]. The method is extended to rectangular pistons as follows. The frequency domain representation of ( ) is obtained from the inverse discrete Fourier transform (DFT) as ∑̂ for , where ̂ After substituting ̂ ̂( ( ( )( ) (2.21) ) is a sampled frequency domain signal with , , ), the result can be expressed as 20 ( ), and . ( ) where ( ∑ ̂( ) ) is the sampling frequency, frequency component, ̂ ( velocity, ( ( is the angular frequency of the th ) is the th frequency component of the transducer surface ) is the th time sample of the transducer surface velocity, and number of time samples. The ) is the ( ) function removes the periodic replicas in (2.21) so that only the original signal ( ) that begins at ( (2.22) ) and ends at is retained. with ( Replacing ( ) ) in (2.22) yields ( ) ∑ ̂( ) ( ( )) (2.23) Comparing (2.23) to the expression for standard Time-Space Decomposition given in Equation (2.18) yields ( ) ( ) The number of terms in (2.18) is ( ̂( ) ) ( (2.24) ) . The signal ( ) is now decomposed into a finite number of unweighted and weighted complex exponential functions ( ) and ( ) which only depend on the spatial term and the temporal term , respectively. Applying (2.24) to (2.19), the transient pressure field generated by a rectangular piston 21 with an arbitrary surface velocity signal can be now calculated with the Fast Nearfield Method combined with Frequency Domain Time-Space Decomposition. Other differences between Frequency Domain Time-Space Decomposition and TimeSpace Decomposition include 1) the functions ( ) and ( ) are complex in Frequency Domain Time-Space Decomposition and are purely real in Time-Space Decomposition, and 2) the number of terms is the number of discrete signal samples in Frequency Domain Time-Space Decomposition while the number of terms is determined by the signal type in Time-Space Decomposition. Since the DFT can be applied to any discretized signal, Frequency Domain Time-Space Decomposition removes the limitation imposed by Time-Space Decomposition and thereby improves the speed of transient pressure field calculations when combined with Fast Nearfield Method. The typical number of superposition operations in standard Time-Space Decomposition is 2, 6, or 8, depending on the signal type. The number of superposition operations in Frequency Domain Time-Space Decomposition is equal to the number of frequency components in DFT and is usually much larger than that in standard Time-Space Decomposition. In addition, all of the frequency components evaluated in Frequency Domain Time-Space Decomposition are complex-valued. In order to improve the performance of Frequency Domain Time-Space Decomposition, a few specific Fourier transform properties are applied. In particular, for real-valued signals, ̂ ( ̂ ( ) ). Thus, only half of the frequency components are needed to represent the original signal. Further improvement can be achieved with spectral clipping. In spectral 22 clipping, a threshold value is defined for the spectrum. Complex-valued frequency components with amplitudes below the threshold are discarded in the calculation, and only components above the threshold are passed to the Frequency Domain TimeSpace Decomposition routine. With spectral clipping, Frequency Domain Time-Space Decomposition is significantly accelerated, and the error increases by only a small amount. 2.2 B-mode imaging simulation In this section, the impulse response simulation model for a pulse-echo system and the convolution model in terms of transmit and receive pressure signals will be explained. Several algorithms designed to improve the performance for the latter model will also be introduced. 2.2.1 Impulse Response Pulse-echo Model For a pulse-echo imaging system, the convolution integral that describes the relationship between the excitation and the received signal are given in Stepanishen [29]. For an pulse-echo input signal ( ), which includes the transducer excitation and the electromechanical impulse response of the transmit and receive transducers that is already known from either established electromechanical properties or from measurements, the pulse-echo model can be represented as [30] ( ) In equation (2.25), transducer, ( ) (2.25) represents the simulated pulse-echo signal received by the is the speed of sound in the medium, and 23 and are the impulse response of the transmit aperture and receive aperture, respectively. Aliasing effects are observed in the impulse response pulse-echo simulation model when the sampling frequency is too low. If ( ) can be expressed as the convolution between two identical signals ( ) [31] such that ( ) ( ) ( ) (2.26) The simulated pulse-echo signal (2.25) can be then expressed as ( ) where ̇( ) ̇( ) is the density of the medium, ⁄( (2.27) ) is the transient pressure at the observation point generated by the transmit aperture when excited by the signal ( ), and is the transient pressure at the same point generated by the receive aperture for an excitation signal ( ). An equivalent expression for the pulse-echo signal terms of and is thus derived in . When the Fast Nearfiled Method is combined with Time-Space Decomposition or Frequency Domain Time-Space Decomposition, pulse-echo signals for B-mode imaging simulations can be calculated faster and with more accuracy than simulations that use the impulse response simulation model. 2.2.2 Algorithms for B-mode imaging simulations In (2.27), the pressure signals and are generated by the transmit and receive apertures, respectively. One approach calculates pressure signals generated from the 24 entire aperture using Fast Nearfield Method. The convolutions between and are then executed with Fast Fourier Transforms using the open source fftw library (http://www.fftw.org/). Figure 2-3 Aperture Movement. The gray transducer elements represent the active aperture. When the aperture is translated to obtain a new A-line, the pressure signals contributed from same transducer elements are calculated repeatedly with the first simulation algorithm. When the aperture is translated to the center of the next A-line, there are repeated calculations for the pressure signals generated by same transducer elements in the next aperture. The contributions from the transducer elements colored in black in Figure 2-3 are calculated twice with this simulation approach when the aperture (colored in gray) moves from Figure 2-3 (a) to Figure 2-3 (b). Considering the linearity of the pulse-echo system, total transmit and receive pressure signals can be written as the superposition 25 of individual pressure signals generated from each transducer element. The convolution between and can be then be represented as ( ) ( ) (2.28) With Equation (2.28), second algorithm that convolves the signals first and then superposes the results has been developed. In this algorithm, the individual pressure signals generated from each transducer element are calculated first. The convolutions between each pair of transmit and receive pressure signals are then evaluated and stored. Since the results of the convolution are the same when transmit and receive apertures are swapped, the resulting waveforms can be saved in an upper-triangular "matrix" structure as (2.29) ( ) When the time delays are digitized for a specific system clock, the RF data can then be calculated by reading these stored waveforms and superposing them with appropriate time delays according to the location of the focal point along each A-line. In Equation (2.28), if the convolutions between each pair of transmit and receive pressure signals are not evaluated and stored and instead individual pressure signals generated from each single element are evaluated and stored, the total aperture 26 pressure signals can then be obtained by superposing individual element pressure signals with appropriate time delays. Convolutions between total transmit and receive aperture pressure signals are then evaluated with FFTs to obtain the RF data. This third algorithm therefore adds the signals first and then performs a single convolution to represent received pulse-echo signal. 2.3 Error Evaluation In calculations of transmitted ultrasound signals in a plane, the normalized error is evaluated with √∑ ( ) ( ) (2.30) √∑ where ( ( ) ) is the pressure amplitude at the spatial coordinates ( a specified instant in time and at a given distance, and ( ) evaluated at ) is the reference pressure calculated by the impulse response method, which is calculated with a high sampling frequency in the same locations. The indices spatial point ( ). and and indicate a particular represent the number of grid points in the and directions, respectively. For the simulation of B-mode imaging, the normalized error for the RF signal is evaluated with 27 ( √∑ ) ( ) (2.31) ( √∑ where ( ) ) represents a single time point on a single A-line, and ( ) indicates the reference pulse-echo signal calculated by the impulse response simulation model at a high sampling frequency. The index the index indicates an individual A-line. indicates a particular time point, and is the number of time points, and is the number of A-lines. The reference RF signals are down-sampled before the comparison with the simulated pulse-echo signals is evaluated. 28 Chapter 3 Simulation Results Results of transient ultrasound propagation and B-mode imaging simulations will be described in this chapter. In section 3.1, the propagation of a transient ultrasound wave is simulated with various models. Simulation results of a transient pressure field generated by a rectangular piston transducer are shown in this section. In section 3.2, B-mode imaging is simulated using the pulse-echo model. Images resulting from the simulations of point target scattering and cyst phantom scattering are shown in this section. All simulations were conducted on a computer with an Intel Core i5 CPU at 2.8GHz, with 8GB RAM. The simulations were run in MATLAB R2010a, 64-bit version. Microsoft Visual C++ 2008 SP1 was used to compile C/C++ MEX routines. 3.1 Transient wave propagation simulations In this section, transient pressure fields will be calculated with different propagation models. The simulation parameters are given in sub-section 3.1.1. The reference pressure field is given in sub-section 3.1.2. The simulation results generated with the impulse response method, Fast Nearfield Method, FNM with Time-Space Decomposition, and the FNM with Frequency Domain Time-Space Decomposition are described in following paragraphs. The calculation times and errors for each propagation model are plotted in figures and compared in the table. 3.1.1 Simulation Parameters The transducer used in the transient ultrasound propagation simulations is a rectangular piston of 1.5mm in width and 2.25mm in height. The center of the transducer is located 29 at the origin of the Cartesian coordinate system. The excitation applied to the transducer surface is a four-cycle hanning weighted sinusoid signal centered at 5MHz with normal velocity amplitude. The transient pressure field is calculated at . The calculated plane is located at from the transducer surface and consists of a 50 x 50 point grid from 0 to the piston width along the x axis and from 0 to the piston height along the y axis. All the simulations are calculated in a lossless medium. The speed of sound . The density of the medium . The transducer and calculation plane are illustrated in Figure 3-1. The rectangular area enclosed by the solid lines in Figure 3-1 indicates the rectangular piston located in the X-Y plane centered at the origin. The rectangular area inside the dashed lines at indicates the calculation plane of the transient pressure field. Figure 3-1 Rectangular transducer piston and simulation coordinates. The rectangular piston is 1.5mm in width and 2.25mm in height. The transient pressure field is calculated in the plane 0.5mm from the transducer surface in the first quadrant. Cartesian coordinates are used for the simulation. 30 3.1.2 Reference pressure field The reference pressure field is calculated with the impulse response at a 10GHz sampling frequency. Figure 3-2 shows the reference pressure field in the plane at . Due to the symmetry of the rectangular transducer, the pressure field is only calculated in the first quadrant. The overall pressure field is symmetric to the coordinate origin of the pressure field shown in Figure 3-2. This quarter of the overall pressure field contains both field regions inside (0 – half width/height) and outside (half width/height – width/height) the transducer’s surface space, and therefore obtains a good representation of the overall transient pressure field. The 10GHz sampling frequency ensures the calculation accuracy of the impulse response used as the reference. 31 Figure 3-2 Reference Pressure Field. The reference pressure field is calculated with the impulse response at a 10GHz sampling frequency. The pressure field is only calculated in the first quadrant due to the symmetry of the transducer. For interpretation of the references to color in this and all other figures the reader is referred to the electronic version of this thesis. 3.1.3 Impulse Response In the simulations using the impulse response model, the pressure field is calculated with sampling frequencies from 100MHz to 4GHz. The calculation error is evaluated by (2.30). Simulation time is recorded for each sampling frequency. Figure 3-3 shows the relationship between calculation time and error for the impulse response model. The axes of the plot are in the log scale. 32 In the impulse response, the pressure field is obtained by a temporal convolution evaluated using the Fast Fourier Transform (FFT) in MATLAB. To optimize the calculation, the signals calculated in the FFT have been zero padded to a power of 2 points. Hence when the sampling frequency increases, the change of the error is relatively stable in some regions and quite sharp in other regions in Figure 3-3. The sudden change in error indicates a change in the number of FFT points used in calculating the convolution. Due to aliasing at the impulse response edge, the calculation error is large when the sampling frequency is low. At a sampling frequency of 100MHz, the impulse response generates a 1.70% calculation error with a computation time of 0.16 seconds. To reach computation accuracy better than 1%, the impulse response requires a sampling frequency of 200MHz. The calculation error is 0.54% when the sampling frequency is 200MHz, and the computation time is 0.39 seconds. For a computation error lower than 0.1%, the impulse response requires a sampling frequency higher than 700MHz. The computation time for a 700MHz sampling frequency is 1.76 seconds. Figure 3-3 illustrates the trend of decreasing error while calculation time increases when the sampling frequency increases in the impulse response calculations. 33 Figure 3-3 Impulse Response Calculation Time vs. Error. At a 100MHz sampling frequency, the calculation error of the impulse response is 1.70%. Increasing the sampling frequency to 200MHz, the calculation error drops to 0.54% with a computation time of 0.39s. A calculation accuracy reaching 0.1% error or less requires a sampling frequency over 700MHz. At a 700MHz sampling frequency, the calculation error is 0.75% and the computation time is 1.76s. When the sampling frequency increases, the calculation time also increases while the calculation error drops. 3.1.4 Time-Space Decomposition In Time-Space Decomposition, the hanning weighted sinusoid signal is divided by the sum of products of 6 pair of sine and cosine functions. The calculation accuracy depends greatly on the accuracy of the Gaussian quadrature. Figure 3-4 shows the relation between calculation error and time as the number of Gaussian abscissas increases. The sampling frequency used in the Time-Space Decomposition calculation 34 is 100MHz. The number of abscissas used for the Gaussian quadrature is from 2 to 60. The axes in the plot are in log scale. In Figure 3-4, when increasing the number of Gaussian abscissas, the Time-Space Decomposition calculation time increases while the calculation error decreases. After 4 Gaussian abscissas, the time-error plot shows a significant drop in calculation error while the increase in calculation time is relatively small. To reach a 1% accuracy, TimeSpace Decomposition requires at least 15 Gaussian abscissas. The corresponding calculation time is 0.05 seconds. To reach 0.1% calculation accuracy, the time-space decomposition requires at least 19 Gaussian abscissas. The corresponding calculation time is 0.07 seconds. After 20 Gaussian abscissas, Time-Space Decomposition gives good accuracy in the transient pressure field calculation. Considering both calculation error and time, 20 abscissas is a good choice for the Time-Space Decomposition method. 35 Figure 3-4 Time-Space Decomposition Calculation Time vs. Error. Calculation accuracy increases when increasing the number of Gaussian abscissas used in the Time-Space Decomposition method. To reach 1% accuracy, Time-Space Decomposition requires 15 abscissas. The computing time for 15 abscissas is 0.05s. 19 abscissas are required to reach 0.1% accuracy. The computing time for 19 abscissas is 0.07s. 3.1.5 Frequency Domain Time-Space Decomposition In Frequency Domain Time-Space Decomposition, both the number of frequency components chosen in spectral clipping and the number of abscissas used in Gaussian quadrature influence the calculation accuracy. The following paragraphs discuss the influence of these two factors in two sections. Section 3.1.5.1 gives the result of spectral clipping and section 3.1.5.2 shows the relation of calculation time and error with and 36 without spectral clipping. The sampling frequency used in Frequency Domain TimeSpace Decomposition calculations is 100MHz. 3.1.5.1 Spectral Clipping In Figure 3-5, four aspects of the four-cycle hanning weighted sinusoid excitation signal are illustrated to explain the spectral clipping. Figure 3-5(a) illustrates the zero-padded original signal. The horizontal axis is the time in microseconds, and the vertical axis is the velocity amplitude. The duration of the excitation signal is 0.8 seconds. When sampled at 100MHz, the number of time points contained in the input signal is 81. To optimize the FFT calculation, the excitation signal has been zero padded to 128 temporal points. Figure 3-5(b) shows the spectrum of the input signal. The horizontal axis is the frequency in MHz, and the vertical axis is the spectrum power in dB. Due to frequency domain symmetry, only half of the spectrum is shown from 0 to 50MHz. The dashed line and the dotted line indicate the power thresholds for 1% and 0.1% clipping at -44dB and -63dB, respectively. Figure 3-5(c) shows the reconstruction error when applying the inverse FFT with different numbers of frequency components. The horizontal axis indicates the number of frequency components used in the reconstruction, and the vertical axis indicates the reconstructed signal error. Since the Fourier transformation property of ̂( ̂ ( ) ) is applied for the real input signal, half of the 128 frequency components can represent all of the information in the spectrum. The reconstruction error is introduced when applying the inverse FFT with fewer frequency components. When clipping is implemented at -44dB, the reconstruction error is below 1%. The 10 frequency components that contain most of the energy in the spectrum are chosen 37 for the reconstruction. To reach reconstruction accuracy higher than 0.1%, the 15 components that contain most of the power are required. The spectral clipping for this accuracy is implemented at -63dB. The dashed and dotted lines in Figure 3-5(c) illustrate the reconstruction error of 1% and 0.1% accuracy with the corresponding number of frequency components, respectively. Figure 3-5(d) shows the reconstructed signal at 0.67% accuracy with the 10 frequency components containing most of the spectral power. The reconstructed signal matches well with the original signal within the duration of the signal with only small oscillations in the zero-padded region. Figure 3-5 Spectral Clipping. (a) The zero-padded original excitation signal. The duration of the original signal is 0.8s with 81 temporal points. The signal is zero-padded to 128 temporal points. (b) Spectrum of the original signal. Half of the frequency domain 38 Figure 3-5 (cont’d) is shown considering the symmetry of the frequency domain. The dashed line and the dotted line indicate the threshold of spectral clipping for 1% and 0.1% accuracy, respectively. (c) Reconstruction error compared to the number of frequency components chosen to reconstruct the signal. The dashed line and dotted line indicate reconstruction errors of 1% and 0.1% accuracy and corresponding number of frequency components, respectively. (d) Reconstructed signal at 1% accuracy. 3.1.5.2 Calculation error – time Figure 3-5 shows the relation of calculation time to error with and without spectral clipping. The calculations use from 2 to 60 Gaussian abscissas. In Figure 3-5, the horizontal axis is the calculation time, and the vertical axis is the calculation error. The horizontal dashed and dotted lines illustrate 1% and 0.1% calculation error, respectively. The solid line indicates the calculation time and error without spectral clipping. All 64 frequency terms are used in this case. The leftmost dashed line indicates the relationship between calculation time and error when spectral clipping is implemented at 1% accuracy, and the dotted line indicates the relationship between calculation time and error when the spectrum clipping is implemented at 0.1% accuracy. The FDTSD method without spectral clipping ultimately reaches an accuracy of 0.001% when the number of Gaussian quadrature abscissas increases to 60. For 1% accuracy, FDTSD without spectral clipping requires 15 Gaussian abscissas and takes 0.44 seconds to finish the calculation. For 0.1% error, 19 Gaussian abscissas are required. The calculation takes 0.54 seconds at this accuracy. For FDTSD with 1% spectral clipping, the simulation result converges to an accuracy of 0.41%. When the number of Gaussian abscissas is large, the calculation error is 39 mainly caused by the 1% spectral clipping. For an error lower than 1%, the 1% spectral clipping FDTSD requires 13 Gaussian abscissas and takes 0.07 seconds to complete the calculation. When spectral clipping is implemented at 1% accuracy, the FDTSD method will not converge to 0.1% calculation accuracy. For FDTSD with 0.1% spectrum clipping, the simulation results converge to 0.096% accuracy in the end. The FDTSD method requires 15 Gaussian abscissas to reach 1% accuracy. The corresponding calculation time is 0.10 seconds. 24 Gaussian abscissas are required to reach 0.1% calculation accuracy which takes 0.15 seconds. Figure 3-6 Frequency Domain Time-Space Decomposition Calculation Time vs. Error. The solid line indicates FDTSD without spectral clipping. To reach 1% accuracy, 40 Figure 3-6 (cont’d) 15 Gaussian abscissas are required and the calculation time is 0.44s. 19 Gaussian abscissas give a calculation accuracy of 0.1%. The corresponding calculation time is 0.54s. The dashed line indicates FDTSD with 1% spectral clipping. 1% calculation accuracy is achieved with 13 Gaussian abscissas with 0.07s computing time. FDTSD with 1% spectral clipping cannot converge to 0.1% accuracy. The dotted line indicates FDTSD with 0.1% spectral clipping. To reach 1% calculation accuracy, the FDTSD with 0.1% spectral clipping requires 15 Gaussian abscissas and 0.10s to complete computation. To reach 0.1% accuracy, the method requires 24 Gaussian abscissas and takes 0.15s. 3.1.6 Comparison Figure 3-6 shows a comparison of the performances of different calculation models. In the figure, the solid line indicate FDTSD without spectral clipping, the dashed line indicates FDTSD with 1% spectral clipping, the dotted line indicates FDTSD with 0.1% spectral clipping, the solid line with circle marks indicates the impulse response, and the solid line with square marks indicates TSD method. In Figure 3-6, calculation time increases to the right along the horizontal axis and calculation error increases upward along the vertical axis. Thus, the time – error lines located to the right of the figure take more calculation time than the lines located to the left of the figure. The horizontal dashed and dotted lines throughout the figure indicate calculation errors of 1% and 0.1%, respectively. Observing the intersection with 1% and 0.1% horizontal lines, the cross-point of the TSD line appears farthest to the left, which indicates that the TSD method takes least calculation time to reach a given accuracy. The next cross-points are the FDTSD method with 1% spectrum clipping, and FDTSD method with 0.1% spectrum clipping. To reach a given accuracy, the FDTSD method executes faster when fewer frequency components are selected in the spectral clipping. For low calculation accuracies, the FDTSD method without spectral clipping can be even slower than the 41 impulse response. When the calculation accuracy increases, the advantage of FDTSD’s calculation speed becomes obvious even in the case of no spectral clipping. At higher accuracies, the impulse response requires high sampling frequencies. The calculation speed of the impulse response with a high accuracy is the slowest of all methods mentioned above. Figure 3-7 Comparison of all propagation models. The horizontal axis is the calculation time, and the vertical axis is the calculation error. From left to right, the solid line with square marks is the TSD method, the dashed line is FDTSD with 1% clipping, the dotted line is FDTSD with 0.1% clipping, the solid line is FDTSD without clipping, and the solid line with circle marks is the impulse response. TSD is located to the far left in the figure, which indicates it is the fastest model in this calculation. The impulse response is located to the far right in the figure when the accuracy is high, which indicates it takes most time to calculate a high accuracy result. 42 Table 3-1 gives an explicit comparison of the calculation time for each propagation model. At 1% accuracy, the Time-Space Decomposition method is about 7x faster than the impulse response. FDTSD without spectral clipping is even slower than the impulse response. FDTSD with spectral clipping can calculate at the same speed as the TSD method. FDTSD with 1% spectral clipping calculates 5x faster than the impulse response, and FDTSD with 0.1% spectral clipping calculates 4x faster than the impulse response. For 0.1% accuracy, the TSD method is 26x faster than the impulse response. FDTSD without spectral clipping calculates 3x faster than the impulse response at this accuracy, and FDTSD with 0.1% spectral clipping calculates 12x faster than the impulse response. FDTSD with 1% spectral clipping cannot reach a 0.1% accuracy in the overall calculation. Table 3-1 Comparison of calculation time Model Impulse Response TSD FDTSD – No Clipping FDTSD – 1% Clipping FDTSD – 0.1% Clipping 3.2 Calculation Time (s) 1% error 0.1% error 0.3924 1.7615 0.0535 0.0677 0.4354 0.5431 0.0743 0.0969 0.1504 B-mode imaging simulations This section will describe B-mode imaging simulations. Sub-section 3.2.1 gives the simulation parameters used in the B-mode imaging simulation. In sub-section 3.2.2, the B-mode imaging is simulated for the point target scattering. Sub-section 3.2.2.1 generates simulations with a single focus point. When no system 43 clock is applied, the simulations are implemented in MATLAB with FOCUS and Field II simulators in 3.2.2.1.A. The algorithms are written in C/C++ and implemented in MATLAB with C/C++ MEX files. The simulation results are shown in 3.2.2.1.B. “Algorithm 1 / ALG1” indicates the first algorithm described in Section 2.2.2. “Algorithm 2 / ALG2” indicates the second algorithm, and “Algorithm 3 / ALG3” indicates the third algorithm. In 3.2.2.1.C, simulations are generated with different sampling frequencies and system clocks to discuss the influence of these two factors on the simulation. Subsection 3.2.2.2 generates simulations with multiple focal zones. In sub-section 3.2.3, the B-mode imaging simulations are generated for the cyst phantom scattering. The simulations are generated with the Field II simulator and with the pulse-echo simulation algorithm 3 (“sum first, convolve second algorithm”). Imaging results and RF data comparisons are shown in this section. 3.2.1 Simulation Parameters In point target scattering simulations, 20 point scatterers are equally arranged from 15mm from the transducer surface to 110mm along z-axis. The distance between scatterers is 5mm. The transducer array used in the point target scattering simulation is composed of 128 rectangular elements. Each rectangular element is 0.5133mm in width and 5mm in height. The kerf between elements is 0.1mm wide. The excitation signal ( ) is obtained by convolving two single cycle sine functions with 3MHz center frequency. In Field II, 2x20 sub-rectangles are applied on each transducer element. The B-mode image is 20mm in width and contains 20 amplitude-lines (A-lines). For each Aline, a sub-aperture of 64 elements is defined centered at the A-line. For simulations with single focal zone, a single focus point is defined at 60mm from aperture surface 44 along each A-aline. For simulations with multiple focal zones, 5 focal zones are defined for each transmit and receive sub-aperture. The focus point of each focal zone locates equally from 30mm to 110mm from transducer surface. The distance between two focus points is 20mm. In Field II, 5 set of time delays are set to sub-apertures independently in beamforming to assure the time delays are identical to those used in the reference. The simulation results are generated in each focal zone and stitched together to compose the Field II results. The simulation times for Field II are recorded for the typical routines generating the multiple-focal-zone results. The B-mode images are generated with 60dB dynamic range. In cyst phantom scattering simulations, 100,000 scatterers are randomly generated within a 50mm x 10mm x 60mm cubic space. The phantom starts at 30mm from the transducer surface. 5 cysts with 0 amplitude are generated at x = 10mm and y = 0 along the z-axis. The centers of cysts are from z = 40mm to 80mm, and with 10mm distance between each. The diameters of 6 cysts are 6mm, 5mm, 4mm, 3mm, and 2mm, respectively. Five highly scattering cysts with 10x amplitude are generated at x = -5mm and y = 0 along the z-axis. The centers of the cysts are from z = 40mm to 80mm, with 10mm distance between each. The diameters of the highly scattering cysts are 2mm, 3mm, 4mm, 5mm, and 6mm. Six point scatterers with 20x amplitudes are generated at x = -15mm and y = 0 along the z-axis. The first point scatterer is located at 40mm from transducer surface. The distance between two point scatterers is 10mm. The transducer array used in cyst phantom simulations is composed of 192 rectangular elements. Transducer elements’ sizes are the same as defined in the point target scattering simulations. The kerf between elements is 0.05mm wide. The excitation signal 45 ( ) is obtained by the convolving two single cycle sine functions with 3.5MHz centering frequency. In Field II, 1x10 sub-rectangles are applied on each transducer element. The B-mode image is 40mm in width and contains 50 amplitude-lines (A-lines). For each Aline, a sub-aperture of 64 elements is defined. A single focus point located at 60mm along each A-line is defined for each aperture. A 40MHz system clock is applied in simulations. The sound speed . The density of the medium . The reference RF data and B-mode image are generated by impulse response at 5GHz sampling frequency. 3.2.2 Point target scattering 3.2.2.1 Single focal zone for transmit and receive apertures 3.2.2.1.A B-mode imaging simulation in FOCUS without system clock Figure 3-8 shows the 5GHz reference impulse response B-mode imaging simulation result. Figure 3-9, Figure 3-10 and Figure 3-11 show comparisons of Field II and FOCUS simulated image results generated at different sampling frequencies. The left columns in these three comparison images are simulated images generated by Field II. The right columns in the comparison images are simulated images generated by FOCUS. In Figure 3-8, 20 point scatterers are shown from 15mm to 110mm in the middle Aline of the simulated image. A single focus point is shown at 60mm. The reference image has a clear view throughout the image due to the high sampling frequency. Only slight blurring appears in the nearfield region. 46 Figure 3-8 Reference B-mode simulated image for single focal zone simulation without system clock. The reference RF data is generated by 5GHz impulse response. 20 point scatterers are shown in the middle of the image from 15mm to 110 mm. The single focal zone is shown at 60mm. 47 Figure 3-9 Field II and FOCUS simulated image comparison at 25MHz sampling frequency. Clutters are obvious in the nearfield region in Field II simulated image, whereas the artifacts are not obvious in FOCUS simulated image. The calculation error for Field II is 18.58% and for FOCUS is 1.47%. Simulation time for Field II is 0.23s and for FOCUS is 1.08s. When the sampling frequency is 25MHz, the simulation error for Field II result is 18.58%, while the simulation error for FOCUS result is 1.47%. In Figure 3-9, clutters 48 appear in the nearfield region in the Field II simulated image, whereas nearfield artifacts are not obvious in the FOCUS simulated image. The calculation time for Field II at 25MHz sampling frequency is 0.23 seconds. For FOCUS, the simulation takes 1.08 seconds to complete computing. Figure 3-10 Field II and FOCUS simulated image comparison at 50MHz sampling frequency. Nearfield artifacts are reduced in the Field II simulated image. The FOCUS simulated image remains high quality. The calculation error for Field II at 50MHz 49 Figure 3-10 (cont’d) sampling frequency is 7.02% and for FOCUS is 0.37%. The calculation time for Field II is 0.26s and for FOCUS is 1.60s. When increasing the sampling frequency to 50MHz, the simulation error for Field II decreases to 7.02%, and the error for FOCUS decreases to 0.37%. Nearfield clutter is reduced in Field II simulated image as shown in Figure 3-10. The FOCUS simulation keeps the high quality in the simulated image. At as low sampling frequency as 50MHz, FOCUS reaches a simulation error lower than 1%. The calculation time for Field II at a 50MHz sampling frequency is 0.26 seconds. The calculation time for FOCUS is 1.60 seconds. 50 Figure 3-11 Field II and FOCUS simulated image comparison at 100MHz sampling frequency. Field II generates a visually clear image at 100MHz and FOCUS keeps the high simulation quality. The calculation error for Field II is 5.07% and for FOCUS is 0.19%. The simulation time for Field II is 0.33s and for FOCUS is 2.62s. At 100MHz sampling frequency, both Field II and FOCUS simulate visually clear result images. The simulation error of Field II is 5.07% and of FOCUS is 0.19%. It 51 takes 0.33 seconds for Field II to finish the simulation, whereas the calculation time for FOCUS is 2.62 seconds at 100MHz sampling frequency. For the FOCUS simulation routine, it only requires 30MHz sampling frequency to reach a simulation error of less than 1%. For Field II, 250MHz sampling frequency and 2x20 spatial sampling are required to reach a simulation error below 1%. The main factors slowing down the FOCUS calculation speed are repeating calculations of transducer elements and the evaluations of FFTs. Thus, several implementing algorithms run in C/C++ routines are developed as following. 3.2.2.1.B Implementing algorithms with 40MHz system clock To take advantage of the system linearity, a 40MHz system clock is applied to the simulation. The time delays for each aperture are digitized and aligned to the system clock. Figure 3-12 shows the reference simulated image generated by the impulse response at 5GHz sampling frequency with 40MHz system clock. Figure 3-13 shows the comparison of simulated images generated by the first algorithm (ALG1), the second algorithm (ALG2) and the third algorithm (ALG3) at 40MHz sampling frequency. In Figure 3-12, artifacts appear in the nearfield region and beyond the focus region due to the digitized time delays. 52 Figure 3-12 Reference B-mode simulated image for single focal zone simulation with 40MHz system clock. The reference RF data are generated by impulse response at 5GHz sampling frequency. Due to the time delays, clutter appears in the nearfield region and beyond the focus region. 53 Figure 3-13 Comparison of implementing algorithms. Three columns of simulated images from left to right are generated by the first algorithm, the second algorithm and the third algorithm, respectively. The time delays are aligned to a system clock of 40MHz. The sampling frequency is 40MHz. The simulation results of all algorithms obtain simulation accuracies of 0.56%. The simulation times are 0.41s, 0.47s and 0.05s for the first algorithm, the second algorithm, and the third algorithm, respectively. Figure 3-13 illustrates the performances of three algorithms at 40MHz sampling frequency. Since all algorithms are based on the same simulation approach and the 54 only difference is the calculation sequence, the simulation error is 0.56% for all three algorithms. For the first algorithm calculating with Fast Nearfield Method with TimeSpace Decomposition, the calculation time is 0.41 seconds. The second algorithm takes 0.47 seconds to complete the simulation, and the third algorithm takes 0.05 seconds to finish the simulation. In the first algorithm, the number of convolutions equals to the number of A-lines multiplied by the number of scatterers. In the second algorithm, the number of convolutions is proportional to the number of elements. In the third algorithm, the number of convolutions is the same as in the first algorithm. Considering the time for evaluating FFTs to calculate convolutions, the second algorithm can be slower than other two algorithms when the number of elements is large. The third algorithm is a good combination of the first algorithm and the second algorithm, and thus speeds up the computation about 8x faster than the first algorithm and the second algorithm. When applying 40MHz sampling frequency and 40MHz digitized time delays to Field II, the simulation takes 0.28 seconds to complete the computation. The simulation error at this low sampling frequency is 15.56% for Field II simulation result, which is much larger than Fast Nearfield Method based approaches. 3.2.2.1.C Influence of sampling frequency and system clock When the digitized time delays are applied in the simulation, both the sampling frequency and the system clock to which the digitized time delays are aligned to will affect the simulation performance. Figure 3-14 and Figure 3-15 illustrate how these two factors influence the simulation. In Figure 3-14, the digitized time delays are aligned to a 40MHz system clock. The simulation results are generated with 55 sampling frequencies of 40MHz, 200MHz and 400MHz and the corresponding simulated images are shown in three columns from left to right in Figure 3-14. In Figure 3-15, the sampling frequency is set to 400MHz. The digitized time delays are aligned to system clocks of 40MHz, 200MHz and 400MHz. The corresponding simulated images are shown in three columns from left to right in Figure 3-15. The reference used in this sub-section is generated by the impulse response at 5GHz sampling frequency without system clock. The simulated image for the reference is shown in Figure 3-8 in section 3.2.2.1.A. When the digitized time delays are aligned to a system clock of 40MHz, the simulation error is 13.77% regardless of the sampling frequency. In Figure 3-14, the cluttering in the nearfield regions in all simulated images reflects the influence of a low system clock. When the system clock is relatively low in the simulation, increasing the sampling frequency cannot generate a simulation result of higher accuracy. In this condition, the system clock is crucial to the calculation accuracy and the quality of simulated image. 56 Figure 3-14 Comparison of simulated images generated with different sampling frequencies. The digitized time delays are aligned to a 40MHz system clock. The three columns of simulated images from left to right are generated with 40MHz, 200MHz and 400MHz sampling frequencies, respectively. The nearfield cluttering indicates the influence of a low system clock. The simulation error is 13.77% regardless of the sampling frequency. When the sampling frequency is 400MHz, the comparison of simulated images generated with different system clocks are shown in Figure 3-15. The left column in 57 Figure 3-15 is the simulated image where the digitized time delays are aligned to a 40MHz system clock. The nearfeild cluttering is caused by the low system clock as described above. When increasing the system clock to 200MHz, the simulation error drops to 2.98%. The middle column in Figure 3-15 shows the simulated image. A visually clear image is generated with this system clock. A system clock of 400MHz with 400MHz sampling frequency can provide a simulation accuracy of 1.44% and create a simulated image as shown in the right column in Figure 3-15. The simulated image obtains a good quality in this condition and is visually comparable to the reference result without system clock. 58 Figure 3-15 Comparison of simulated images generated with different system clocks. The sampling frequencies are 400MHz for all simulations. The three columns of simulated images from left to right are generated with 40MHz, 200MHz and 400MHz system clocks, respectively. The nearfield cluttering is obvious when the system clock is low as 40MHz. When the system clock increases to 200MHz and 400MHz, the artifacts visually disappear. The simulation errors are 13.77%, 2.98% and 1.44% for 40MHz, 200MHz and 400MHz system clocks, respectively. 59 Table 3-2 gives a conclusion of the influences of sampling frequencies and system clocks. The system clock to which the digitized time delays are aligned is a crucial factor affecting the simulation quality. A system clock of 200MHz can provide a visually acceptable quality in the simulated image. Table 3-2 Sampling frequency and system clock effect 40MHz system clock Sampling frequency Error System Clock Error fs = 40MHz fs = 200MHz 13.77% 13.77% 400MHz sampling frequency Clk = 40MHz Clk = 200MHz 13.77% 2.89% fs = 400MHz 13.77% Clk = 400MHz 1.44% 3.2.2.2 Multiple focal zones for transmit and receive aperture Figure 3-16 shows the reference B-mode simulated image for 5 transmit and receive focal zones generated by impulse response at 5GHz with 40MHz system clock. Figure 3-17 shows the comparison of simulated images generated by Field II and the third algorithm at 40MHz sampling frequency with 40MHz system clock. Figure 3-18 shows the comparison of simulated images generated by Field II and the third algorithm at 40MHz sampling frequency with 40MHz system clock. 60 Figure 3-16 Reference B-mode simulated image for 5 transmit and receive zones. The reference result is generated by the impulse response at 5GHz sampling frequency. The digitized time delays are aligned to a 40MHz system clock.Nearfield cluttering is cause by the low system clock. The multiple focal zones improves the simulation quality in focal regions. In Figure 3-16, the cluttering in the nearfield region and beyond the focal regions is caused by the system clock. Due to the multiple focal zones, the artifacts are 61 reduced in the focal regions. Since the threshold of two focal zones are set in the middle of two focus points, some edge affects are shown at the thresholds (40mm, 60mm, etc.). At a sampling frequency of 40MHz, the Field II simulation obtains a simulation error of 16.64%. The calculation time for the typical Field II multi-focal-zone routine is 0.31 seconds. For the third algorithm, the simulation error is 0.45% at 40MHz sampling frequency. The simulation completes in 0.08 seconds. Increasing the sampling frequency to 200MHz, the typical Field II routine reaches an accuracy of 1.96%. The computation time for the typical Field II routine at 200MHz sampling frequency is 0.51 seconds. The third algorithm obtains a simulation accuracy of 0.19% at 200MHz sampling frequency. The computation time for the third algorithm is 0.24 seconds. At a low sampling frequency as 40MHz, the third algorithm already reaches a high simulation accuracy of 0.45%, which is even better than the simulation accuracy obtained by Field II at 200MHz sampling frequency. The computing speed of the third algorithm is faster than the typical Field II simulation routine as well. 62 Figure 3-17 Comparison of Field II and third algorithm at 40MHz sampling frequency with 40MHz system clock. The left column is the simulated image generated by Field II at and the right column is the simulated image generated by the third algorithm (ALG3). Field II reaches a simulation accuracy of 16.64% in 0.31s, whereas algorithm 3 reaches a simulation accuracy of 0.45% in 0.08s. 63 Figure 3-18 Comparison of Field II and third algorithm at 200MHz sampling frequency with 40MHz sampling frequency. The left column is the simulated image generated by Field II at and the right column is the simulated image generated by the third algorithm (ALG3). Field II reaches a simulation accuracy of 1.96% in 0.51s, whereas algorithm 3 reaches a simulation accuracy of 0.19% in 0.24s. 64 Table 3-3 Point target scattering B-mode imaging simulation comparison Simulation Type Single Focal zone 5 Focal Zones Sampling Frequency No system clock 25MHz 50MHz 100MHz 40MHz System Clock 40MHz 40MHz System Clock 40MHz 200MHz Approach Error Time Error Time Field II FOCUS 18.58% 0.23s 1.47% 1.08s 7.02% 0.26s 0.37% 1.60s 5.07% 0.33s 0.19% 2.62s 15.56% 0.28s ALG1 ALG2 ALG3 0.56% 0.41s 0.56% 0.47s 0.56% 0.05s Field II ALG3 16.64% 0.31s 0.45% 0.08s 1.96% 0.51s 0.19% 0.24s Table 3-3 summarizes the performances of point target scattering simulations in section 3.2.2. For the simulations with a single focal zone and without system clock, at same sampling frequencies, the FOCUS simulations can generate results 13x 27x more accurate than Field II but with longer computation times. When the time delays are digitized and aligned to a 40MHz system clock, the third algorithm (ALG3) using the Fast Nearfield Method with Time-Space Decomposition generates simulation results with the same accuracy but 8x faster than the first algorithm (ALG1). In the simulations with multiple focal zones, the third algorithm (ALG3) generates simulation results that are 37x more accurate than Field II at 40MHz sampling frequency with 40MHz system clock. The computation time for ALG3 is 4x faster than Field II. At 200MHz sampling frequency with 40MHz system clock, the ALG3 generates 10x more accurate results than Field II and the computing time is 65 2x faster than Field II. The simulation results generated by the third algorithm (ALG3) at 40MHz sampling frequency are even more accurate than Field II results at 200MHz sampling frequency. 3.2.3 Cyst Phantom scattering Figure 3-19 and Figure 3-20 show the simulated images of cyst phantom generated by Field II and the third algorithm (ALG3), respectively. Zero amplitude cysts, highly scattering regions, and point scatterers can be seen in the simulated images. Since only a single focal zone is applied in these simulations, the cysts beyond the focal region can be obscured or deformed. From the discussions of simulations of point target scattering, a sampling frequency of 200MHz is appropriate for Field II simulations to obtain an acceptable accuracy. To complete a simulation with 100,000 random scatterers, Field II requires more than 1 hour to finish the computation. For the third algorithm (ALG3), a sampling frequency of 40MHz is enough to obtain accurate results. The third algorithm only needs approximately 7 minutes to finish the simulation. The RF results generated by Field II and the third algorithm are compared in Figure 3-21. 66 Figure 3-19 Cyst phantom simulated image generated by Field II. The simulation results are generated at 200MHz sampling frequency with a 40MHz system clock. Five zero amplitude cysts can be seen along x = 10mm, five highly scattering regions can be seen along x = -5mm, and 6 point scatters can be seen along x = -15mm. The image regions within in the focal zone are clearer and more accurate than the regions beyond the focal zone. The simulation time for Field II at 200MHz sampling frequency is more than 1 hour. 67 Figure 3-20 Cyst phantom simulated image generated by the third algorithm (ALG3). The simulation results are generated at 40MHz sampling frequency with a 40MHz system clock. Five zero amplitude cysts can be seen along x = 10mm, five highly scattering regions can be seen along x = -5mm, and 6 point scatters can be seen along x = -15mm. The image regions within in the focal zone are clearer and more accurate than the regions beyond the focal zone. The simulation time for ALG3 at 40MHz sampling frequency is around 7 minutes. In Figure 3-21, RF data along the 7th A-line (x = -15.2mm) generated by Field II and the third algorithm (ALG3) are compared. Five point scatterers are located near this A-line. In Figure 3-21, the horizontal axis is the axial distance along z-axis, and the vertical axis is the RF data amplitude normalized by the maximal signal value. Zooming in 58mm to 62mm along the axial distance, the signal of the 3rd scatterer located at 60mm is 68 distinct from the background. As shown in Figure 3-21, RF results generated by the third algorithm at 40MHz sampling frequency (marked as circles in the figure) match well with the results generated by Field II (the solid line in the figure) at 200MHz sampling frequency. Figure 3-21 Cyst Phantom RF data comparison. The horizontal axis is the axial distance along z-axis, the vertical axis is the signal amplitude normalized by the maximal signal amplitude. The solid lines indicates the RF results generated by Field II at 200MHz sampling frequency. The circle marks are the RF results generated by the third algorithm (ALG3) at 40MHz sampling frequency. RF results generated by ALG3 at 20MHz sampling frequency match well with the results generated by Field II at 200MHz sampling frequency. 69 Chapter 4 4.1 Discussions Transient wave propagation models Different propagation models are used for different purposes in ultrasound simulations. The Rayleigh-Sommerfeld integral has analytical solutions and is often used in continuous wave calculations. Due to its slow convergence rate in transient pressure calculations, this method is less commonly used in the transient pressure field simulations. For transient pressure field calculations, the impulse response method, Fast Nearfield Method, Time-Space Decomposition method and Frequency Domain Time-Space Decomposition method have been calculated and discussed in this thesis. 4.1.1 Impulse Response The impulse response is one of the most popular calculation models for transient pressure propagation. Derived from the Green’s function, the impulse response method simplifies the calculation by one convolution between the time derivative of the velocity and the spatial impulse response. Closed form spatial impulse response functions for different geometries can be obtained by integrals given in Equation (2.7) and have been derived in several papers. Due to the discontinuities at the edges in the impulse response function, a high sampling frequency is required in the calculation to obtain an accurate result. At high sampling frequencies, the impulse response can generate highly accurate results, yet the computing time will be long. The impulse response result generated at a high sampling frequency is often used as a reference. 70 4.1.2 Fast Nearfield Method Fast Nearfield Method reduces the number of integrals evaluated in the RayleighSommerfeld integral method, isolates the repeated calculations, combines integrals sharing same integrand and common upper or lower limits, and therefore speeds up the computation. Furthermore, since no edge discontinuity occurs in the Fast Nearfield Method, the aliasing problem that occurs in the impulse response model at low sampling frequencies has been avoided. As a consequence, Fast Nearfield Method can reach a high accuracy at a low sampling frequency. Considering the velocity containing both temporal and spatial terms in the transient pressure field calculations, Fast Nearfield Method is often used with Time-Space Decomposition or Frequency Domain Time-Space Decomposition. 4.1.3 Time-Space Decomposition Time-Space Decomposition optimizes the Fast Nearfield Method by dividing the velocity term into the superposition of temporal and spatial functions. The expression for temporal and spatial functions for specific signal types can be found in Table 2-1. As shown in section 3.1.6, the calculation speed of Fast Nearfield Method with Time-Space Decomposition is the fastest in all computing methods discussed in this thesis. The limitation for Time-Space Decomposition is that the expressions for temporal and spatial functions can be obtained only for several signal types. Thus the Frequency Domain Time-Space Decomposition has been developed to take advantage of the decomposition. 71 4.1.4 Frequency Domain Time-Space Decomposition Frequency Domain Time-Space Decomposition applies the decomposition to the velocity signal on arbitrary signal types. With spectral clipping, the trade-off between computing time and calculation accuracy can be controlled by selecting appropriate number of frequency components. As shown in section 3.1.6, the Frequency Domain Time-Space Decomposition method can generate accurate result at a comparable speed to the TSD method when 1% spectral clipping is applied. The calculation time increases when accuracy of the spectral clipping increases. The FDTSD is an ideal calculation model when the excitation signal cannot be decomposed by TSD method (such as the Gaussian-modulated sinusoidal pulse). The computing time of FDTSD method will be much smaller than the impulse response method when the spectral clipping is applied appropriately. 4.2 B-mode imaging simulations The impulse response model for a pulse-echo system can be simplified as equation (2.25). Simulated results generated by impulse response at a high sampling frequency are used for reference in this thesis. By expressing the excitation signal in a convolutional form, the impulse response model can be written as the convolution between transmit and receive pressure signals. The RF data can be then calculated with efficient transient pressure method as discussed in the previous section. Two types of B-mode imaging simulations are generated in this thesis: the point target scattering and the cyst phantom scattering. 72 4.2.1 Point target scattering 4.2.1.1 Single Focal Zone In simulation software FOCUS, the transient pressure signals are calculated by Fast Nearfield Method with Time-Space Decomposition using existing routines. Implementing the convolutional model given in (2.27) directly in FOCUS, the simulation reaches higher accuracy compared to Field II at the same sampling frequency. To reach same accuracies, FOCUS requires much lower sampling frequencies than the impulse response based software Field II. Several implementing algorithms are developed to improve the performance of the FNM based imaging simulation routine. With digitized time delays that are aligned to a specific system clock, the first algorithm implements the convolutional model directly from equation (2.27). Repeated calculations are executed for individual pressure signals generated from same transducer elements. The second algorithm calculates the convolution between each pair of transmit and receive individual pressure signals and stores the convolutional waveform. The RF data are generated by these stored waveforms with appropriate time delays. The third algorithm calculates and stores each individual pressure signal. Instead of convolving between each pair of individual signals, the third algorithm applies time delays to the stored individual signals and superposes them to generate aperture pressure signals. The convolutions are then executed between transmit and receive aperture signals. From comparison in Table 3-3, the second algorithm takes the longest computing time, and the third algorithm simulates the fastest. Considering the number of convolutions evaluated in each algorithm, the first algorithm and the third algorithm evaluate the number of convolutions equal to the product of the number of A-lines 73 multiplying the number of scatterers ( ). The second algorithm, in the other hand, evaluates the convolutions proportional to the square of the number of transducer elements times the number of scatterers ( ) ( ). As a consequence, when the number of elements in the transducer is huge, the computing time for this algorithm will be long. The influence of sampling frequencies and system clocks have been considered. When the digitized time delays are aligned to a low system clock, increasing the sampling frequency will not improve the simulation performance. A system clock of 200MHz provides acceptable simulation quality. For all implementing algorithms, the requirement of a sampling frequency that is higher than the system clock needs to be satisfied. 4.2.1.2 Multiple focal zones B-mode imaging simulation algorithms can be applied for simulations with multiple focal zones. Thresholds for different focal zones are set at the middle point between two focus points. The multiple focal zones improve the performance in the focal regions. Artifacts exist at thresholds due to the edge problem. The third algorithm provides an implementing method to generate accurate results with a low sampling frequency requirement for B-mode imaging simulations with multiple focal zones. 4.2.2 Cyst Phantom scattering Simulations for complicated images with a large number of scatterers can be generated by B-mode imaging simulation algorithms discussed in the previous section. Since the third algorithm has the best performance in the point target scattering simulations, we implement this simulation algorithm in the 100,000-scatterer cyst phantom scattering 74 imaging simulation. At a sampling frequency of 40MHz, the third algorithm generates RF results which are well matching to the results generated by Field II at 200MHz. The third algorithm takes only 7 minutes to finish the 100,000-scatterer imaging simulation, which is 9x faster than the Field II simulation. For complicated imaging simulations, the third algorithm shows its potential to be a fast and accurate simulation approach. 75 Chapter 5 Conclusions In calculations of transient pressure field, the impulse response method generates high accuracy results with high sampling frequencies. The Fast Nearfield Method avoids the aliasing problem in the impulse response method and therefore reaches the same accuracies with lower sampling frequency requirements. The Time-Space Decomposition further speeds up the Fast Nearfield Method by dividing temporal and spatial terms and combing repeated calculations. To break the restriction of signal type in Time-Space Decomposition, Frequency Domain Time-Space Decomposition has been developed. Spectral Clipping helps FDTSD to control the trade-off between calculation accuracy and computing time. For B-mode imaging simulations, the convolutional expression in equation (2.27) introduces an approach applying FNM with TSD/FDTSD in B-mode imaging simulations. The FNM based imaging simulation approach requires lower sampling frequencies to reach the same simulation accuracies. The third algorithm generates accurate imaging results with low sampling frequency requirements and high computing speeds. The system clock to which the digitized time delays are aligned will affect the simulation performance. The implementing algorithms can be applied for imaging simulations with single focal zone or with multiple focal zones. The speed advantage of the third algorithm stands out when applied to simulations for complicated images. 76 Chapter 6 Future Work In the B-mode imaging simulations with multiple focal zones, better windows for separating different focal zones can be developed to avoid artifacts at the zones’ edges. The implementing routine for the cyst phantom imaging simulation with multiple focal zones can be developed with the third algorithm to obtain better simulation results. Apodization on sub-apertures can further improve the simulation quality. With the development of convolutional waveform expressions and the FDTSD method, the limitation of the excitation signal type can be reduced. Good interpolation approaches should also be considered to improve the simulation performance when the system clock is limiting the simulation quality. 77 REFERENCES 78 REFERENCES [1]. C. I. Zanelli, C. W. Hennige, N. T. Sanghvi, “Design and characterization of a 10 cm annular array transducer for high intensity focused ultrasound (HIFU) applications”, Ultrasonics Symposium, 1994. Proceedings., 1994 IEEE , 3: 1887 – 1890. (1994) [2]. R. J. McGough, T. V. Samulski, and J. F. Kelly. “An efficient grid sectoring method for calculations of the nearfield pressure generated by a circular piston.” J. Acoust. Soc. Am., 115(5):1942–1954 (2004). [3]. J. Hoffelner, M. Bechtold, B. Granz, R. Lerch, K. Newerla, "Self-focusing HIFU source for large therapy volumes," Ultrasonics Symposium, 1998. Proceedings., 1998 IEEE , vol.2, no., pp.1563,1566 vol.2, 1998 [4]. J. Georgii, C. V. Dresky, S. Meier, D. Demedts, C. Schumann, and T. Preusser. "Focused Ultrasound-Efficient GPU Simulation Methods for Therapy Planning." In Proc. Workshop on Virtual Reality Interaction and Physical Simulation VRIPHYS, Lyon, France, pp. 119-128. (2011). [5]. K. B. Ocheltree and L. A. Frizzell, Sound field calculation for rectangular sources, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 36: 242–248 (1989). [6]. J. C. Lockwood and J. G. Willette, High-speed method for computing the exact solution for the pressure variations in the nearfield of a baffled piston. J. Acoust. Soc. Am., 53(3): 735-741 (1973). [7]. F. Oberhettinger, On transient solutions of the “baffled piston” problem. Journal of Research of the National Bureau of Standards - B. Mathematics and Mathematical Physics, 65B(1) : 1- 6 (1961) [8]. J. A. Jensen, Ultrasound fields from triangular apertures. J. Acoust. Soc. Am., 100(4): 2049-2056 (1996). [9]. A. Penttinen and M. Luukkala, The impulse response and pressure nearfield of a curved ultrasonic radiator. J. Phys. D: Appl. Phys, 9(10):1547–1557 (1976). [10]. P. Faure, D. Cathignol, J. Y. Chapelon, and V. L. Newhouse. On the pressure field of a transducer in the form of a curved strip. J. Acoust. Soc. Am., 95(2):628– 637 (1994). [11]. J. A. Ketterling. Acoustic field of a wedge-shaped section of a spherical cap transducer. J. Acoust. Soc. Am., 114:3065–3075 (2003). 79 [12]. X. Zeng, R. J. McGough, Evaluation of the angular spectrum approach for simulations of near-field pressures, J. Acoust. Soc. Am. Volume 123, Issue 1, pp. 68-76 (2008) [13]. X. Zeng, R. J. McGough, Optimal simulations of ultrasonic fields produced by large thermal therapy arrays using the angular spectrum approach, J. Acoust. Soc. Am. Volume 125, Issue 5, pp. 2967-2977 (2009) [14]. U. Vyas, D. Christensen, Ultrasound Beam Propagation using the Hybrid Angular Spectrum Method, Engineering in Medicine and Biology Society, 2008. EMBS 2008. 30th Annual International Conference of the IEEE , vol., no., pp.2526,2529, 20-25 Aug. 2008 [15]. Robert J. McGough, Rapid calculations of time-harmonic nearfield pressures produced by rectangular pistons, J. Acoust. Soc. Am., 115(5): 1934-1941 (2004). [16]. J. Huang, R. G. Holt, R. O. Cleveland, and R. A. Roy, Experimental validation of a tractable numerical model for focused ultrasound heating in flow-through tissue phantoms. J. Acoust. Soc. Am. 116, 2451 (2004) [17]. G. Kossoff, Basic Physics and Imaging Characteristics of Ultrasound, World J. Surg., 24: 134-142. (2000) [18]. Peter R. Stepanishen, Transient Radiation from Pistons in an Infinite Planar Baffle, J. Acoust. Soc. Am., 49(5B): 1629-1638 (1971). [19]. L. G. Ullate and J. L. San Emeterio, A new algorithm to calculate the transient near-field of ultrasonic phased-arrays. IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 39(6): 745–753 (2001). [20]. M. Arditi, F. S. Foster, and J. W. Hunt. Transient fields of concave annular arrays. Ultrason. Imaging, 3(1):37–61 (1981). [21]. J. A. Jensen, N. B. Svendsen, "Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers," Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on , vol.39, no.2, pp.262,267, March 1992 [22]. J. A. Jensen, "Simulating arbitrary-geometry ultrasound transducers using triangles," Ultrasonics Symposium, 1996. Proceedings., 1996 IEEE , vol.2, no., pp.885,888 vol.2, 3-6 Nov 1996 [23]. J. D’hooge, J. Nuyts, B. Bijnens, B. De Man, P. Suetens, J. Thoen, M.-C. Herregods, and F. Van de Werf. The calculation of the transient near and far field of a baffled piston using low sampling frequencies. J. Acoust. Soc. Am., 102(1):78–86 (1997). 80 [24]. J. A. Jensen, "A new calculation procedure for spatial impulse responses in ultrasound", J. Acoust. Soc. Am. 105, 3266 (1999) [25]. J. F. Kelly and R. J. McGough, A Time-Space Decomposition Method for Calculating the Nearfield Pressure Generated by a Pulsed Circular Piston, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 53(6): 1150-1159 (2006). [26]. D. Chen and R. J. McGough. A 2D fast near-eld method for calculating nearfield pressures generated by apodized rectangular pistons. J. Acoust. Soc. Am., 124: 1526-1537 (2008). [27]. E. J. Alles, Y. Zhu, K. W. A. van Dongen, R. J. McGough, Rapid Transient Pressure Field Computations in the Nearfield of Circular Transducers Using Frequency-Doma [28]. V. Chan and A. Perlas, “Basics of Ultrasound Imaging”, Atalas of UltrasoundGuided Procedures in Interventional Pain Management, pp 13-19, Springer Science + Business Media, LLC. (2011) [29]. Peter R. Stepanishen, Pulsed transmit/receive response of ultrasonic piezoelectric transducers, J. Acoust. Soc. Am., 69(6): 1815-1827 (1981). [30]. J. A. Jensen, A Model for the Propagation and Scattering of Ultrasound in Tissue, J. Acoust. Soc. Am., 89(1): 182-191 (1991). [31]. Y. Zhu, R. J. McGough, and T. L. Szabo, "A comparison of ultrasound image simulations with FOCUS and Field II," Proc. 2012 IEEE Ultrasonics Symposium, Dresden, Germany, 2012. [32]. J. A. Jensen, Field: A Program for Simulating Ultrasound Systems, Med. Biol. Eng. Comp., 10th Nordic-Baltic Conference on Biomedical Imaging, Suppl.1, Part 1, 34: 351-353 (1996). [33]. G. E. Tupholme, "Generation of acoustic pulses by baffled plane pistons", Mathematika, 16:209-224, 1969. 81