- Email: [email protected]
Optimization of point spread function in ultrasound arrays
Ultrasonics 44 (2006) 159–165 www.elsevier.com/locate/ultras
Optimization of point spread function in ultrasound arrays M. Sakhaei a, A. Mahloojifar a
a,*
, A. Malek
b
Department of Biomedical Engineering, Tarbiat Modarres University, P.O. Box 14115-143 Tehran, Iran b Department of Mathematics, Tarbiat Modarres University, P.O. Box 14115-175 Tehran, Iran Received 13 April 2005; received in revised form 26 September 2005; accepted 13 October 2005 Available online 15 November 2005
Abstract The design of aperture weightings in ultrasound arrays is a multi-objective optimization problem, involving parameters such as delays, aperture size, focal depth, operating frequency and beam properties. Besides, apodization causes the SNR in array output to be decreased. We introduce an analytic expression of the lateral point spread function and a model for SNR as nonlinear functions of weights, based on which, a new aperture design method is established, resulting in an optimal set of weights. These weights provide a point spread function having the predetermined peak sidelobe level, while the SNR in array output is optimized. Optimization results from a linear array with M = 128 elements equally spaced at one wavelength, center frequency f0 = 3.5 MHz and 50% relative bandwidth, have shown that decreasing the peak sidelobe level, decreases the SNR. Therefore, an array designer can select a proper set of weights according to its application, using a SNR curve versus to the peak sidelobe level. In addition, the method can maintain the same beam properties over a long range with low variations in the SNR. Simulation results have shown only 1 dB variations in the SNR for depths from 20 mm to 120 mm, which is a longer range and better SNR performance than the conventional methods. Ó 2005 Elsevier B.V. All rights reserved. Keywords: Linear array; Apodization; Point spread function; Sensitivity; Resolution
1. Introduction In diagnostic ultrasound, contrast resolution and lateral resolution are two main measures of image quality. They affect the ability to detect deep lying, low contrast lesions in the body, such as tumors [1]. In order to increase lateral resolution, it is essential that the beam width in the focal point is no wider than a couple of wavelengths of the imaging frequency. For the contrast resolution, the beam intensity should diminish quickly beyond the focal point. One of the simple and common ways to control the beam properties and image quality is apodization, i.e. weighting of the array elements [2]. To find the optimal weights, various methods have been proposed. Some methods deal with narrowband beam pattern [1–4] and the others with wideband
*
Corresponding author. E-mail addresses: [email protected] (M. Sakhaei), [email protected] (A. Mahloojifar). 0041-624X/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2005.10.002
beam pattern [5–10]. Besides, there are some researches concerned with finding the optimal weights to obtain the desired beam properties at various depths [11–13]. In addition, different criteria have been proposed for the objective function. The maximization of the ratio of main lobe energy to sidelobe energy [2], the minimization of the peak sidelobe level [3,6,5], the minimization of the sidelobe energy [6,11] are some introduced objective functions. Another objective function is minimization of the error between the obtained and the desired function over sidelobe region [4,7,9], or over main lobe and sidelobe regions [10,12,13]. The constant delay and the constant excitation amplitude weight in each channel are usually used for transmission of the wave. However, in receive mode, dynamic receive focusing may be used, by which the channel delay continually changes with time [14]. This provides the best focusing at each depth. Additionally, the aperture length may also change dynamically so that the f-number stays constant. The f-number is the ratio of the range being
160
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
interrogated to the aperture size. This method, named dynamic receive aperture, uses fewer elements for shallower depths, and therefore may not be an optimal method for image enhancement, particularly in the presence of noise. Moreover, in unequally spaced arrays, the method cannot be implemented as well as in equally spaced arrays. One of the undesired effects of apodization is the reduction of the signal to noise ratio (SNR) in the array output. In many practical situations, the SNR is reduced in imaging of deep lying tissue and the increased reduction of SNR by apodization may decrease depth of field. This becomes critical in two-dimensional arrays which have smaller elements and lower sensitivity in comparison with the elements of one-dimensional arrays [15]. In this paper, we propose a new method to calculate the optimal weights in the process of designing the desired beam properties. The method is based on a new analytic expression for the one-way lateral PSF and a model analyzing the effect of apodization on SNR. To design the desired PSF, an optimization problem to control the beam characteristics, while minimizing the corruptive effects of apodization is proposed and solved. Our approach is similar to the technique proposed by Holm and Elgeton [3], but they optimized continuous wave beam pattern and not the wideband PSF in near field. Furthermore, there are differences between our method and that described by Ranganathan and Walker [12]. Our proposed method is able to control the sidelobe level more precisely. In addition, neither of them applies any constraint on weights to optimize the SNR. The present paper is organized as follows. In Section 2, an analytic expression for the one-way lateral PSF is established, then the mathematical expression for the SNR as a function of the weights is described. The new mathematical optimization problem is given in Section 3. In Section 4, the performance of the proposed method is evaluated through some examples. Conclusions and further discussions are presented in the last section.
z F(0,zf)
-d/2
P(x,zf)
xk d/2
x
Fig. 1. Geometry used for calculating lateral point spread function.
element directivity effects, the output signal of the array can be approximated by gðx; tÞ ¼
M X
wk sðt T ðx; xk ÞÞ
ð1Þ
k¼1
where wk is the weight of the kth element, s(t) is the pulse wave emanated from the source, T(x, xk) is the summation of path delay (Td) and beamforming delay (Tb). For the homogenous medium, Td and Tb are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 Td ¼ ðxk xÞ þ z2 ð2-1Þ c sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f 0 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 d2 þ z2f x2k þ z2f A Tb ¼ @ ð2-2Þ c 4 where c is the wave propagation velocity in the medium, d is the aperture length and the phase reference is assumed to be on one side of the array. The envelope of the array output can be calculated by env ½gðx; tÞ ¼ jgðx; tÞ þ j^ gðx; tÞj
ð3Þ
where g^ðx; tÞ is the time domain Hilbert transform of g(x, t). Since Hilbert transform is a linear operator, using Eq. (1) in Eq. (3), it can be rewritten as X M env ½gðx; tÞ ¼ wk ðsðt T ðx; xk ÞÞ þ j^sðt T ðx; xk ÞÞÞ k¼1 ð4Þ
2. Theory 2.1. Lateral point spread function An ultrasound imaging system is usually considered as a linear system [16]. Therefore, it is possible to characterize beam properties by the point spread function or spatial impulse response. A mathematical expression of the receive lateral PSF is here described and that of the transmit lateral PSF can be established in a similar way. Consider an array having M = 2N elements in the xz plane where the origin is in the center of the array. The x-axis lies on the line connecting elements, while z is vertical to the array plane, as shown in Fig. 1. Assume that the array is focused on a point F at the depth of zF on the z-axis (main axis) and a point source lies at point P(x, zF) in the neighborhood of point F(0, zF). Thus, by ignoring the attenuation of the medium, beam spreading and the
We know that s(t) is a pass band signal that may be modeled as a carrier with frequency f0 amplitude modulated with a low pass signal. Assuming m(t) as the low pass signal, Eq. (4) will be approximated as X M env ½gðx; tÞ ¼ wk mðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞ k¼1 ð5Þ jxt
jxt
in which the relation Hilbert {m(t) e } = m(t)j e is employed. We use Eq. (6) to define the lateral point spread function in point x psf ðxÞ ¼ Max jgðx; tÞj t
ð6Þ
which is similar to that considered in [5–10], for calculating wideband beam pattern. This expression is convenient in
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
which is the time that the wave reaches to the phase reference. Therefore, the PSF may be represented as in Eq. (8) X M ð8Þ wk mðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞ psf ðxÞ ¼ k¼1 This equation is in analytic form and explicitly represents lateral PSF as a nonlinear function of weights and the pulse wave. In addition, it reduces the computational burden needed for calculating PSF, compared to Eq. (6). For symmetric arrays and symmetric weights, where p(x, t) = p(x, t), it is possible to reduce the computational burden more gðx; tÞ ¼
N X
0
Normalized Magnitude (dB)
-5 -10 -15 -20 -25 -30 -35 -40 -45 -50 -20
-15
-10
-5
0
5
10
15
20
Lateral Position (mm) (a) 0
-10
Normalized Magnitude (dB)
cases when the optimization problem can be converted to linear programming. In a linear array, the lateral point spread function around the focal point spreads on a line parallel to the array [16]. Therefore, the maximum of env [g(x, t)] at x is the value of array output corresponding to point F(0, zF) assuming a point source on P(x, zF) (see Fig. 1). In the other words, if a wave travels from point F at time t = 0, the value of point spread function occurs at the time it reaches to the phase reference element. Since the wave is assumed to be originated from point P at t = 0, therefore the psf (x) can be obtained by the value of env [g(x, t)] at time t = t*, where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 d 2 þ jxj t ¼ zf þ ð7Þ c 2
161
-20
-30
-40
-50
-60
wk ðsðt T ðx; kÞÞ þ sðt T ðx; kÞÞÞ
ð9Þ
k¼1
and psf (x) may be written as X N w ðmðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞ psf ðxÞ ¼ k¼1 k þ mðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞÞ
sðtÞ ¼ expððpBf 0 tÞ =aÞ cosð2pf0 tÞ
-15
-10
-50
0
5
10
15
20
Lateral Position (mm) (b)
Fig. 2. Comparison of calculated (solid line) and simulated (dotted line) PSF for two set of weightings: (a) uniform and (b) hamming.
ð10Þ
To show that t* considered in Eq. (7) is properly selected and to evaluate the Eq. (10), the ultrasound field produced by a linear array is simulated by FIELDII [17] and its PSF is computed according to the definition in (6). Then the results are compared to that obtained by Eq. (10). The array has M = 128 elements equally spaced at one wavelength (corresponding to the center frequency), element width L = 0.2 mm and center frequency f0 = 3.5 MHz and is focused on point F(0,60 mm). The pulse wave is assumed as follows: 2
-20
ð11Þ
where a = 1.2 Ln10 and B, the 6 dB relative bandwidth, is set to 0.5. Fig. 2 shows the simulated and calculated PSFs for two sets of weightings, uniform weightings and Hamming
weightings. In two cases, the PSF is computed over 201 equally spaced lateral points between 20 mm and 20 mm. As Fig. 2 shows, Eqs. (7) and (10) are in good agreement with the simulation results. 2.2. Signal to noise ratio analysis Let us assume that a point source is at the focal point of an array. Thus, the array output signal can be obtained as X yðtÞ ¼ wk rk ðtÞ ð12Þ k
where rk(t) is the signal at kth channel after applying beamforming delay and can be expressed as rk ðtÞ ¼ sk ðtÞ þ gk ðtÞ
ð13Þ
where sk(t) is the main signal and gk(t) is the transducer noise. As the source is assumed on the focal point, sk(t)Õs are equal and y(t) can be written as
162
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
X
yðtÞ ¼
! wk sðtÞ þ
k
X
wk gk ðtÞ
ð14Þ
spread function for x P xc is not higher than de times of the PSF at the central point.
k
If gk(t)Õs are independent and have the same power of r2n , then signal to noise ratio in output will be equal to P 2 wk ð15Þ SNR ¼ Pk 2 SNR0 k wk where SNR0 is the signal to noise in each channel after beamforming delay. However, from CauchyÕs inequality, we know that P 2 wk Pk 2 6 M ð16Þ k wk where M is the number of elements and equality is valid only when the wkÕs are the same. This means that in the uniform apodization case, the SNR value at array output is M times of the SNR in a channel and apodization decreases the SNR value. Hence, to increase the SNR value, M should be increased and the weights distribution should be close to uniform. It should be noted that the ratio of SNR in array output to SNR in each channel, before beamforming delay is referred to array gain, and in presence of spatially uncorrelated noise, is called white noise gain. Both are frequency dependent in wideband cases [18]. 3. Optimization In this section, we express and solve the problem of optimally designing the point spread function, considering the corruptive effects of apodization on SNR. Our main objective is to find the weights so that the PSF has desired properties and at the same time, the SNR is maximized. As stated above, apodization can reduce the magnitude of PSF at sidelobes i.e. far points from the main axis and therefore enhance the image contrast resolution. Considering poor sensitivity of elements, the design of optimal weights may be expressed as the following: Min
N X
w2k
k¼1
Subject to N X wk ¼ 1 2
ð17Þ
k¼1
psfðxÞ 6 d e
for xc 6 x 6 xend
where xend is the maximum image width, xc is a control point that can determine the main lobe width. The array and its weighting are also assumed symmetric. According to Eq. (15), minimizing P the sum of the squared weights subN ject to the constraint 2 k¼1 wk ¼ 1 means maximizing the SNR. Therefore, this approach provides the weights to maximize signal to noise such that the magnitude of point
Fig. 3. Optimization results of an array described in the text for different values of de: (a) point spread function for de = 20, 30, . . . , 90 dB, (b) weights curve and (c) signal to noise ratio.
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
The optimization problem has a second order objective function and second order constraints (using Eq. (10) for psf(x)), which may be solved by nonlinear programming methods [19]. We solve the problem using fmincon function in optimization toolbox of MATLAB6.5. The minimum possible value for de can be obtained from the optimization problem in (18), so that the problem stated in (17) has no feasible solution for de lower than that obtained from the problem (18) Min e Subject to N X 2 wk ¼ 1
ð18Þ
k¼1
psfðxÞ 6 e for xc 6 x 6 xend This problem is similar to that considered by Holm and Elgeton [3] for beam pattern function but not for the point spread function. This form does not consider the SNR problem caused by weighting and is therefore convenient in cases where the sensitivity is not a restricting problem. 4. Simulation results In this section the results of the optimization problem (17) for the array described in Section 2.1 are presented. The problem is solved for xc = 2 mm, xend = 20 mm and some values for de. The results obtained for focal depth zF = 60 mm are shown in Fig. 3. The figure demonstrates that smaller sidelobe level can be obtained by increasing the dynamic range of weights and consequently decreasing the SNR. According to Fig. 3, the method may find the weights such that by lowering the de from 20 dB to 80 dB, SNR is only decreased by 5.5 dB. More decreas-
163
ing the de, forces SNR to a very low value, which is not desirable. Therefore, it seems there is a threshold for the maximum of sidelobe level, decreasing from which, needs to loss the SNR greatly. To investigate the changes of the point spread function in different depths, the optimization problem is solved for depths zf = 20, 40, 80, 100, 120 and 140 mm, assuming xc = 2 mm, xend = 20 mm and de = 40 dB. The results are shown in Fig. 4. It is obvious that the weights can be calculated at each depth such that the peak sidelobe level of PSF remains unchanged when increasing the depth from 20 mm to 120 mm while the SNR value varies only about 1 dB. Fig. 4 also shows that extending the properties of the PSF to depth of zf = 140 mm, causes the SNR to decrease rapidly to zero. This means, as for de, there is a threshold level for the range of depths, in which the PSF can have the same maximum sidelobe level. Extending the range beyond the threshold diminishes the SNR quickly. In addition, the method can be employed to implement the optimal dynamic apodization, as illustrated in Fig. 4a. To compare the efficacy of the method to that of conventional dynamic apodization, the PSFÕs obtained by dynamic apodization with dynamic focusing is depicted in Fig. 5a. To this end, we used Hamming window at the depth of zf = 60 mm and allowed the apodization profile to grow as a function of depth and also assuming infinitely large aperture. Then the central portion of the profile is considered to implement the dynamic apodization in 128 receive elements. For depths less than 60 mm, windowing is applied on less receive elements to maintain constant F-number [20]. Figs. 4 and 5b show the SNR and the magnitude of sidelobe levels (at xc = 2 mm) obtained by dynamic apodization at different depths. Comparing them to those obtained by our method shows that our method
Fig. 4. Optimization results in different focal depths for xc = 2 mm and de = 40 dB: (a) point spread function for zf = 20, 40, 60, 80, 100, 120 and 140 mm and (b) the obtained SNR versus zf (solid line) and that obtained by dynamic apodization (dashed line).
164
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
Fig. 5. Conventional dynamic apodization: (a) point spread function for zf = 20, 40, 60, 80, 100, 120 and 140 mm and (b) the magnitude of PSF at xc = 2 mm versus zf (solid line) and that obtained by our method (dashed line).
can more effectively maintain the same peak sidelobe level in depths from 20 mm to 120 mm while the SNR values are similar. 5. Conclusions Apodization of array elements is an effective approach to control the sidelobe level of the point spread function and also to enhance contrast resolution. On the other hand, due to the poor sensitivity of the elements, apodization reduces signal to noise ratio. In this paper, a new method to find the optimal weights for enhancement of the contrast resolution and also to control the signal to noise ratio is proposed. The method provides the weights to produce a PSF having desired properties while maximizing the SNR. It employs an analytic expression of the one-way wideband lateral PSF and a model representing the SNR as functions of weights. The optimal weights are calculated through solving a mathematical optimization problem with the goal to obtain a PSF whose sidelobe magnitude is lower than a predetermined value and the objective to minimize the sum of squared weights as a measure of SNR. In no apodization case, the SNR reaches its maximum possible value, which is equal to the number of elements, and constraining the PSF, causes the SNR to be decreased. Therefore, as simulation results show, there is a trade-off between the SNR and constraints, such that an array designer can select a proper set of weights according to its applications using a SNR curve versus constraints (see Fig. 3c or Fig. 4b). The method has clear advantages in 2D arrays, where the sensitivity of array elements is usually lower than in 1D array. Therefore, weighting in 2D array should be care-
fully applied and it is necessary to control the reduction of SNR by weighting, as the method does. An important application of the method is the optimal dynamic apodization. The optimization method could be used to find the optimal weights in dynamic apodization and as simulation results have shown, it is more effective to maintain the beam characteristics in deep layers in comparison with the conventional dynamic apodization. Moreover, it may be appropriately applied to unequally spaced arrays. It should be noted that the noise discussed in this paper is the internal noise of beamformer system, consisting thermal noise of transducer and electronic system and also random errors from the implementation of apodization. In addition, it may be approximately extended to the phase errors emanating from implementation of delays and from non-homogeneity in the medium. Moreover, due to the linearity of the system in the range of the gain variations, the noise power is assumed to be independent of the system gain. However the optimization problem tries to reduce the sum of squared weights and consequently p reduces ffiffiffiffiffiffiffiffiffiffiffi P 2 the magnitude of the largest weight (i.e. jwmax j 6 wk ), resulting in reduction of the range of gain variations. Finally, the optimization method presented here, may also be applied to two-way point spread function to find optimal receive weights, which have the potential of enhancement the depth of field, because it tries to maintain the same beam properties over depth while maximizing the SNR, both of which, are important to enhance the depth of field. In addition, using FIR filter in each channel may create better performance than simple amplitude weighting, while decreasing received and internal noise. Therefore,
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
extending the problem to deal with FIR filters is considered as the future work. References [1] P.C. Li, S.W. Flax, E.S. Ebbini, M. OÕDonnel, Blocked element compensation in phased array imaging, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 40 (1993) 283–292. [2] C. Boni, M. Richard, S. Barbarossa, Optimal configuration and weighting of nonuniform arrays according to a maximum ISLR criterion, in: IEEE Int. Conf. Acoust., Speech, Sign. Proc., vol. V, 1994, pp. 57–160. [3] S. Holm, B. Elgetun, G. Dahl, Properties of the beam pattern of weight and layout-optimized sparse arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 44 (5) (1997) 983–991. [4] A. Trucco, Thinning and weighting of large planar arrays by simulated annealing, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 46 (2) (1999) 347–355. [5] G. Cardone, C. Cincotti, P. Gori, M. Pappalardo, Optimization of wide band linear arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 48 (4) (2001) 943–952. [6] G. Cardone, C. Cincotti, P. Gori, M. Pappalardo, Design of wideband arrays for low side-lobe level beam patterns by simulated annealing, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 49 (8) (2002) 1050–1059. [7] A. Trucco, Synthesizing wide-band sparse arrays by simulated annealing, in: Proc. MTS/IEEE Conference and Exhibition on OCEANS, vol. 2, 2001, pp. 989–994. [8] J.L. Schwartz, B.D. Steinberg, Ultrasparse, ultrawideband arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 45 (2) (1998) 376– 393.
165
[9] A. Trucco, S. Curletto, Flattening the side-lobes of wide-band beam patterns, IEEE J. Ocean. Eng. 28 (4) (2003) 760–762. [10] S. Curletto, A. Trucco, On the shaping of the main lobe in wide-band arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 52 (4) (2005) 619–630. [11] O. Keitmann-Curdes, B. Brendel, C. Marg, H. Ermert, Optimization of apodization based on the sidelobe pressure energy in simulated ultrasound fields, in: Proc. IEEE Int. Ultrason. Symp., 2002, pp. 1677–1680. [12] K. Ranganathan, W.F. Walker, A novel beamformer design method for medical ultrasound. Part I: theory, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 50 (Jan.) (2003) 15–24. [13] K. Ranganathan, W.F. Walker, A novel beamformer design method for medical ultrasound. Part II: simulation, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 50 (Jan.) (2003) 25–39. [14] K.E. Thomenius, Evolution of ultrasound beamformers, in: Proc. IEEE Int. Ultrason. Symp., 1996, pp. 1615–1622. [15] S.W. Smith, R.E. Davidsen, C.D. Emery, R.L. Goldberg, E.D. Light, Update on 2D Array Transducers for Medical Ultrasound, in: Proc. IEEE Int. Ultrasonics Symposium, 1995, pp. 1273–1278. [16] R.Z. Zemp, C.K. Abbey, M.F. Insana, Linear system models for ultrasonic imaging, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 50 (Jun.) (2003) 642–654. [17] J.A. Jensen, Field: a program for simulating ultrasound systems, Med. Biol. Eng. Comp. 10th Nordic-Baltic Conference on Biomedical Imaging 4 (Suppl. 1) (1996) 351–353, Part 1. [18] H. Cox, R.M. Zeskind, T. Kooij, Practical supergain, IEEE Trans. Acoust., Speech, Signal Process. ASSP-34 (June) (1986) 393–398. [19] S.S. Rao, Optimization, Theory and Applications, Wiley Eastern Limited, 1984. [20] C. Daft, W. Engeler, Windowing of wide-band ultrasound transducers, in: Proc. IEEE Ultrason. Symp., 1996, pp. 1541–1544.
Optimization of point spread function in ultrasound arrays M. Sakhaei a, A. Mahloojifar a
a,*
, A. Malek
b
Department of Biomedical Engineering, Tarbiat Modarres University, P.O. Box 14115-143 Tehran, Iran b Department of Mathematics, Tarbiat Modarres University, P.O. Box 14115-175 Tehran, Iran Received 13 April 2005; received in revised form 26 September 2005; accepted 13 October 2005 Available online 15 November 2005
Abstract The design of aperture weightings in ultrasound arrays is a multi-objective optimization problem, involving parameters such as delays, aperture size, focal depth, operating frequency and beam properties. Besides, apodization causes the SNR in array output to be decreased. We introduce an analytic expression of the lateral point spread function and a model for SNR as nonlinear functions of weights, based on which, a new aperture design method is established, resulting in an optimal set of weights. These weights provide a point spread function having the predetermined peak sidelobe level, while the SNR in array output is optimized. Optimization results from a linear array with M = 128 elements equally spaced at one wavelength, center frequency f0 = 3.5 MHz and 50% relative bandwidth, have shown that decreasing the peak sidelobe level, decreases the SNR. Therefore, an array designer can select a proper set of weights according to its application, using a SNR curve versus to the peak sidelobe level. In addition, the method can maintain the same beam properties over a long range with low variations in the SNR. Simulation results have shown only 1 dB variations in the SNR for depths from 20 mm to 120 mm, which is a longer range and better SNR performance than the conventional methods. Ó 2005 Elsevier B.V. All rights reserved. Keywords: Linear array; Apodization; Point spread function; Sensitivity; Resolution
1. Introduction In diagnostic ultrasound, contrast resolution and lateral resolution are two main measures of image quality. They affect the ability to detect deep lying, low contrast lesions in the body, such as tumors [1]. In order to increase lateral resolution, it is essential that the beam width in the focal point is no wider than a couple of wavelengths of the imaging frequency. For the contrast resolution, the beam intensity should diminish quickly beyond the focal point. One of the simple and common ways to control the beam properties and image quality is apodization, i.e. weighting of the array elements [2]. To find the optimal weights, various methods have been proposed. Some methods deal with narrowband beam pattern [1–4] and the others with wideband
*
Corresponding author. E-mail addresses: [email protected] (M. Sakhaei), [email protected] (A. Mahloojifar). 0041-624X/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2005.10.002
beam pattern [5–10]. Besides, there are some researches concerned with finding the optimal weights to obtain the desired beam properties at various depths [11–13]. In addition, different criteria have been proposed for the objective function. The maximization of the ratio of main lobe energy to sidelobe energy [2], the minimization of the peak sidelobe level [3,6,5], the minimization of the sidelobe energy [6,11] are some introduced objective functions. Another objective function is minimization of the error between the obtained and the desired function over sidelobe region [4,7,9], or over main lobe and sidelobe regions [10,12,13]. The constant delay and the constant excitation amplitude weight in each channel are usually used for transmission of the wave. However, in receive mode, dynamic receive focusing may be used, by which the channel delay continually changes with time [14]. This provides the best focusing at each depth. Additionally, the aperture length may also change dynamically so that the f-number stays constant. The f-number is the ratio of the range being
160
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
interrogated to the aperture size. This method, named dynamic receive aperture, uses fewer elements for shallower depths, and therefore may not be an optimal method for image enhancement, particularly in the presence of noise. Moreover, in unequally spaced arrays, the method cannot be implemented as well as in equally spaced arrays. One of the undesired effects of apodization is the reduction of the signal to noise ratio (SNR) in the array output. In many practical situations, the SNR is reduced in imaging of deep lying tissue and the increased reduction of SNR by apodization may decrease depth of field. This becomes critical in two-dimensional arrays which have smaller elements and lower sensitivity in comparison with the elements of one-dimensional arrays [15]. In this paper, we propose a new method to calculate the optimal weights in the process of designing the desired beam properties. The method is based on a new analytic expression for the one-way lateral PSF and a model analyzing the effect of apodization on SNR. To design the desired PSF, an optimization problem to control the beam characteristics, while minimizing the corruptive effects of apodization is proposed and solved. Our approach is similar to the technique proposed by Holm and Elgeton [3], but they optimized continuous wave beam pattern and not the wideband PSF in near field. Furthermore, there are differences between our method and that described by Ranganathan and Walker [12]. Our proposed method is able to control the sidelobe level more precisely. In addition, neither of them applies any constraint on weights to optimize the SNR. The present paper is organized as follows. In Section 2, an analytic expression for the one-way lateral PSF is established, then the mathematical expression for the SNR as a function of the weights is described. The new mathematical optimization problem is given in Section 3. In Section 4, the performance of the proposed method is evaluated through some examples. Conclusions and further discussions are presented in the last section.
z F(0,zf)
-d/2
P(x,zf)
xk d/2
x
Fig. 1. Geometry used for calculating lateral point spread function.
element directivity effects, the output signal of the array can be approximated by gðx; tÞ ¼
M X
wk sðt T ðx; xk ÞÞ
ð1Þ
k¼1
where wk is the weight of the kth element, s(t) is the pulse wave emanated from the source, T(x, xk) is the summation of path delay (Td) and beamforming delay (Tb). For the homogenous medium, Td and Tb are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 Td ¼ ðxk xÞ þ z2 ð2-1Þ c sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f 0 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 d2 þ z2f x2k þ z2f A Tb ¼ @ ð2-2Þ c 4 where c is the wave propagation velocity in the medium, d is the aperture length and the phase reference is assumed to be on one side of the array. The envelope of the array output can be calculated by env ½gðx; tÞ ¼ jgðx; tÞ þ j^ gðx; tÞj
ð3Þ
where g^ðx; tÞ is the time domain Hilbert transform of g(x, t). Since Hilbert transform is a linear operator, using Eq. (1) in Eq. (3), it can be rewritten as X M env ½gðx; tÞ ¼ wk ðsðt T ðx; xk ÞÞ þ j^sðt T ðx; xk ÞÞÞ k¼1 ð4Þ
2. Theory 2.1. Lateral point spread function An ultrasound imaging system is usually considered as a linear system [16]. Therefore, it is possible to characterize beam properties by the point spread function or spatial impulse response. A mathematical expression of the receive lateral PSF is here described and that of the transmit lateral PSF can be established in a similar way. Consider an array having M = 2N elements in the xz plane where the origin is in the center of the array. The x-axis lies on the line connecting elements, while z is vertical to the array plane, as shown in Fig. 1. Assume that the array is focused on a point F at the depth of zF on the z-axis (main axis) and a point source lies at point P(x, zF) in the neighborhood of point F(0, zF). Thus, by ignoring the attenuation of the medium, beam spreading and the
We know that s(t) is a pass band signal that may be modeled as a carrier with frequency f0 amplitude modulated with a low pass signal. Assuming m(t) as the low pass signal, Eq. (4) will be approximated as X M env ½gðx; tÞ ¼ wk mðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞ k¼1 ð5Þ jxt
jxt
in which the relation Hilbert {m(t) e } = m(t)j e is employed. We use Eq. (6) to define the lateral point spread function in point x psf ðxÞ ¼ Max jgðx; tÞj t
ð6Þ
which is similar to that considered in [5–10], for calculating wideband beam pattern. This expression is convenient in
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
which is the time that the wave reaches to the phase reference. Therefore, the PSF may be represented as in Eq. (8) X M ð8Þ wk mðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞ psf ðxÞ ¼ k¼1 This equation is in analytic form and explicitly represents lateral PSF as a nonlinear function of weights and the pulse wave. In addition, it reduces the computational burden needed for calculating PSF, compared to Eq. (6). For symmetric arrays and symmetric weights, where p(x, t) = p(x, t), it is possible to reduce the computational burden more gðx; tÞ ¼
N X
0
Normalized Magnitude (dB)
-5 -10 -15 -20 -25 -30 -35 -40 -45 -50 -20
-15
-10
-5
0
5
10
15
20
Lateral Position (mm) (a) 0
-10
Normalized Magnitude (dB)
cases when the optimization problem can be converted to linear programming. In a linear array, the lateral point spread function around the focal point spreads on a line parallel to the array [16]. Therefore, the maximum of env [g(x, t)] at x is the value of array output corresponding to point F(0, zF) assuming a point source on P(x, zF) (see Fig. 1). In the other words, if a wave travels from point F at time t = 0, the value of point spread function occurs at the time it reaches to the phase reference element. Since the wave is assumed to be originated from point P at t = 0, therefore the psf (x) can be obtained by the value of env [g(x, t)] at time t = t*, where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 d 2 þ jxj t ¼ zf þ ð7Þ c 2
161
-20
-30
-40
-50
-60
wk ðsðt T ðx; kÞÞ þ sðt T ðx; kÞÞÞ
ð9Þ
k¼1
and psf (x) may be written as X N w ðmðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞ psf ðxÞ ¼ k¼1 k þ mðt T ðx; xk ÞÞ expðj2pf0 T ðx; xk ÞÞÞ
sðtÞ ¼ expððpBf 0 tÞ =aÞ cosð2pf0 tÞ
-15
-10
-50
0
5
10
15
20
Lateral Position (mm) (b)
Fig. 2. Comparison of calculated (solid line) and simulated (dotted line) PSF for two set of weightings: (a) uniform and (b) hamming.
ð10Þ
To show that t* considered in Eq. (7) is properly selected and to evaluate the Eq. (10), the ultrasound field produced by a linear array is simulated by FIELDII [17] and its PSF is computed according to the definition in (6). Then the results are compared to that obtained by Eq. (10). The array has M = 128 elements equally spaced at one wavelength (corresponding to the center frequency), element width L = 0.2 mm and center frequency f0 = 3.5 MHz and is focused on point F(0,60 mm). The pulse wave is assumed as follows: 2
-20
ð11Þ
where a = 1.2 Ln10 and B, the 6 dB relative bandwidth, is set to 0.5. Fig. 2 shows the simulated and calculated PSFs for two sets of weightings, uniform weightings and Hamming
weightings. In two cases, the PSF is computed over 201 equally spaced lateral points between 20 mm and 20 mm. As Fig. 2 shows, Eqs. (7) and (10) are in good agreement with the simulation results. 2.2. Signal to noise ratio analysis Let us assume that a point source is at the focal point of an array. Thus, the array output signal can be obtained as X yðtÞ ¼ wk rk ðtÞ ð12Þ k
where rk(t) is the signal at kth channel after applying beamforming delay and can be expressed as rk ðtÞ ¼ sk ðtÞ þ gk ðtÞ
ð13Þ
where sk(t) is the main signal and gk(t) is the transducer noise. As the source is assumed on the focal point, sk(t)Õs are equal and y(t) can be written as
162
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
X
yðtÞ ¼
! wk sðtÞ þ
k
X
wk gk ðtÞ
ð14Þ
spread function for x P xc is not higher than de times of the PSF at the central point.
k
If gk(t)Õs are independent and have the same power of r2n , then signal to noise ratio in output will be equal to P 2 wk ð15Þ SNR ¼ Pk 2 SNR0 k wk where SNR0 is the signal to noise in each channel after beamforming delay. However, from CauchyÕs inequality, we know that P 2 wk Pk 2 6 M ð16Þ k wk where M is the number of elements and equality is valid only when the wkÕs are the same. This means that in the uniform apodization case, the SNR value at array output is M times of the SNR in a channel and apodization decreases the SNR value. Hence, to increase the SNR value, M should be increased and the weights distribution should be close to uniform. It should be noted that the ratio of SNR in array output to SNR in each channel, before beamforming delay is referred to array gain, and in presence of spatially uncorrelated noise, is called white noise gain. Both are frequency dependent in wideband cases [18]. 3. Optimization In this section, we express and solve the problem of optimally designing the point spread function, considering the corruptive effects of apodization on SNR. Our main objective is to find the weights so that the PSF has desired properties and at the same time, the SNR is maximized. As stated above, apodization can reduce the magnitude of PSF at sidelobes i.e. far points from the main axis and therefore enhance the image contrast resolution. Considering poor sensitivity of elements, the design of optimal weights may be expressed as the following: Min
N X
w2k
k¼1
Subject to N X wk ¼ 1 2
ð17Þ
k¼1
psfðxÞ 6 d e
for xc 6 x 6 xend
where xend is the maximum image width, xc is a control point that can determine the main lobe width. The array and its weighting are also assumed symmetric. According to Eq. (15), minimizing P the sum of the squared weights subN ject to the constraint 2 k¼1 wk ¼ 1 means maximizing the SNR. Therefore, this approach provides the weights to maximize signal to noise such that the magnitude of point
Fig. 3. Optimization results of an array described in the text for different values of de: (a) point spread function for de = 20, 30, . . . , 90 dB, (b) weights curve and (c) signal to noise ratio.
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
The optimization problem has a second order objective function and second order constraints (using Eq. (10) for psf(x)), which may be solved by nonlinear programming methods [19]. We solve the problem using fmincon function in optimization toolbox of MATLAB6.5. The minimum possible value for de can be obtained from the optimization problem in (18), so that the problem stated in (17) has no feasible solution for de lower than that obtained from the problem (18) Min e Subject to N X 2 wk ¼ 1
ð18Þ
k¼1
psfðxÞ 6 e for xc 6 x 6 xend This problem is similar to that considered by Holm and Elgeton [3] for beam pattern function but not for the point spread function. This form does not consider the SNR problem caused by weighting and is therefore convenient in cases where the sensitivity is not a restricting problem. 4. Simulation results In this section the results of the optimization problem (17) for the array described in Section 2.1 are presented. The problem is solved for xc = 2 mm, xend = 20 mm and some values for de. The results obtained for focal depth zF = 60 mm are shown in Fig. 3. The figure demonstrates that smaller sidelobe level can be obtained by increasing the dynamic range of weights and consequently decreasing the SNR. According to Fig. 3, the method may find the weights such that by lowering the de from 20 dB to 80 dB, SNR is only decreased by 5.5 dB. More decreas-
163
ing the de, forces SNR to a very low value, which is not desirable. Therefore, it seems there is a threshold for the maximum of sidelobe level, decreasing from which, needs to loss the SNR greatly. To investigate the changes of the point spread function in different depths, the optimization problem is solved for depths zf = 20, 40, 80, 100, 120 and 140 mm, assuming xc = 2 mm, xend = 20 mm and de = 40 dB. The results are shown in Fig. 4. It is obvious that the weights can be calculated at each depth such that the peak sidelobe level of PSF remains unchanged when increasing the depth from 20 mm to 120 mm while the SNR value varies only about 1 dB. Fig. 4 also shows that extending the properties of the PSF to depth of zf = 140 mm, causes the SNR to decrease rapidly to zero. This means, as for de, there is a threshold level for the range of depths, in which the PSF can have the same maximum sidelobe level. Extending the range beyond the threshold diminishes the SNR quickly. In addition, the method can be employed to implement the optimal dynamic apodization, as illustrated in Fig. 4a. To compare the efficacy of the method to that of conventional dynamic apodization, the PSFÕs obtained by dynamic apodization with dynamic focusing is depicted in Fig. 5a. To this end, we used Hamming window at the depth of zf = 60 mm and allowed the apodization profile to grow as a function of depth and also assuming infinitely large aperture. Then the central portion of the profile is considered to implement the dynamic apodization in 128 receive elements. For depths less than 60 mm, windowing is applied on less receive elements to maintain constant F-number [20]. Figs. 4 and 5b show the SNR and the magnitude of sidelobe levels (at xc = 2 mm) obtained by dynamic apodization at different depths. Comparing them to those obtained by our method shows that our method
Fig. 4. Optimization results in different focal depths for xc = 2 mm and de = 40 dB: (a) point spread function for zf = 20, 40, 60, 80, 100, 120 and 140 mm and (b) the obtained SNR versus zf (solid line) and that obtained by dynamic apodization (dashed line).
164
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
Fig. 5. Conventional dynamic apodization: (a) point spread function for zf = 20, 40, 60, 80, 100, 120 and 140 mm and (b) the magnitude of PSF at xc = 2 mm versus zf (solid line) and that obtained by our method (dashed line).
can more effectively maintain the same peak sidelobe level in depths from 20 mm to 120 mm while the SNR values are similar. 5. Conclusions Apodization of array elements is an effective approach to control the sidelobe level of the point spread function and also to enhance contrast resolution. On the other hand, due to the poor sensitivity of the elements, apodization reduces signal to noise ratio. In this paper, a new method to find the optimal weights for enhancement of the contrast resolution and also to control the signal to noise ratio is proposed. The method provides the weights to produce a PSF having desired properties while maximizing the SNR. It employs an analytic expression of the one-way wideband lateral PSF and a model representing the SNR as functions of weights. The optimal weights are calculated through solving a mathematical optimization problem with the goal to obtain a PSF whose sidelobe magnitude is lower than a predetermined value and the objective to minimize the sum of squared weights as a measure of SNR. In no apodization case, the SNR reaches its maximum possible value, which is equal to the number of elements, and constraining the PSF, causes the SNR to be decreased. Therefore, as simulation results show, there is a trade-off between the SNR and constraints, such that an array designer can select a proper set of weights according to its applications using a SNR curve versus constraints (see Fig. 3c or Fig. 4b). The method has clear advantages in 2D arrays, where the sensitivity of array elements is usually lower than in 1D array. Therefore, weighting in 2D array should be care-
fully applied and it is necessary to control the reduction of SNR by weighting, as the method does. An important application of the method is the optimal dynamic apodization. The optimization method could be used to find the optimal weights in dynamic apodization and as simulation results have shown, it is more effective to maintain the beam characteristics in deep layers in comparison with the conventional dynamic apodization. Moreover, it may be appropriately applied to unequally spaced arrays. It should be noted that the noise discussed in this paper is the internal noise of beamformer system, consisting thermal noise of transducer and electronic system and also random errors from the implementation of apodization. In addition, it may be approximately extended to the phase errors emanating from implementation of delays and from non-homogeneity in the medium. Moreover, due to the linearity of the system in the range of the gain variations, the noise power is assumed to be independent of the system gain. However the optimization problem tries to reduce the sum of squared weights and consequently p reduces ffiffiffiffiffiffiffiffiffiffiffi P 2 the magnitude of the largest weight (i.e. jwmax j 6 wk ), resulting in reduction of the range of gain variations. Finally, the optimization method presented here, may also be applied to two-way point spread function to find optimal receive weights, which have the potential of enhancement the depth of field, because it tries to maintain the same beam properties over depth while maximizing the SNR, both of which, are important to enhance the depth of field. In addition, using FIR filter in each channel may create better performance than simple amplitude weighting, while decreasing received and internal noise. Therefore,
M. Sakhaei et al. / Ultrasonics 44 (2006) 159–165
extending the problem to deal with FIR filters is considered as the future work. References [1] P.C. Li, S.W. Flax, E.S. Ebbini, M. OÕDonnel, Blocked element compensation in phased array imaging, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 40 (1993) 283–292. [2] C. Boni, M. Richard, S. Barbarossa, Optimal configuration and weighting of nonuniform arrays according to a maximum ISLR criterion, in: IEEE Int. Conf. Acoust., Speech, Sign. Proc., vol. V, 1994, pp. 57–160. [3] S. Holm, B. Elgetun, G. Dahl, Properties of the beam pattern of weight and layout-optimized sparse arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 44 (5) (1997) 983–991. [4] A. Trucco, Thinning and weighting of large planar arrays by simulated annealing, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 46 (2) (1999) 347–355. [5] G. Cardone, C. Cincotti, P. Gori, M. Pappalardo, Optimization of wide band linear arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 48 (4) (2001) 943–952. [6] G. Cardone, C. Cincotti, P. Gori, M. Pappalardo, Design of wideband arrays for low side-lobe level beam patterns by simulated annealing, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 49 (8) (2002) 1050–1059. [7] A. Trucco, Synthesizing wide-band sparse arrays by simulated annealing, in: Proc. MTS/IEEE Conference and Exhibition on OCEANS, vol. 2, 2001, pp. 989–994. [8] J.L. Schwartz, B.D. Steinberg, Ultrasparse, ultrawideband arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 45 (2) (1998) 376– 393.
165
[9] A. Trucco, S. Curletto, Flattening the side-lobes of wide-band beam patterns, IEEE J. Ocean. Eng. 28 (4) (2003) 760–762. [10] S. Curletto, A. Trucco, On the shaping of the main lobe in wide-band arrays, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 52 (4) (2005) 619–630. [11] O. Keitmann-Curdes, B. Brendel, C. Marg, H. Ermert, Optimization of apodization based on the sidelobe pressure energy in simulated ultrasound fields, in: Proc. IEEE Int. Ultrason. Symp., 2002, pp. 1677–1680. [12] K. Ranganathan, W.F. Walker, A novel beamformer design method for medical ultrasound. Part I: theory, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 50 (Jan.) (2003) 15–24. [13] K. Ranganathan, W.F. Walker, A novel beamformer design method for medical ultrasound. Part II: simulation, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 50 (Jan.) (2003) 25–39. [14] K.E. Thomenius, Evolution of ultrasound beamformers, in: Proc. IEEE Int. Ultrason. Symp., 1996, pp. 1615–1622. [15] S.W. Smith, R.E. Davidsen, C.D. Emery, R.L. Goldberg, E.D. Light, Update on 2D Array Transducers for Medical Ultrasound, in: Proc. IEEE Int. Ultrasonics Symposium, 1995, pp. 1273–1278. [16] R.Z. Zemp, C.K. Abbey, M.F. Insana, Linear system models for ultrasonic imaging, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 50 (Jun.) (2003) 642–654. [17] J.A. Jensen, Field: a program for simulating ultrasound systems, Med. Biol. Eng. Comp. 10th Nordic-Baltic Conference on Biomedical Imaging 4 (Suppl. 1) (1996) 351–353, Part 1. [18] H. Cox, R.M. Zeskind, T. Kooij, Practical supergain, IEEE Trans. Acoust., Speech, Signal Process. ASSP-34 (June) (1986) 393–398. [19] S.S. Rao, Optimization, Theory and Applications, Wiley Eastern Limited, 1984. [20] C. Daft, W. Engeler, Windowing of wide-band ultrasound transducers, in: Proc. IEEE Ultrason. Symp., 1996, pp. 1541–1544.