PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Ultrason Ferroelectr Freq Control. Author manuscript; available in PMC 2010 November 8.
Published in final edited form as:
IEEE Trans Ultrason Ferroelectr Freq Control. 2010 June; 57(6): 1248–1262.
PMCID: PMC2975328
NIHMSID: NIHMS214914

Estimation of Scattering Object Characteristics for Image Reconstruction Using a Nonzero Background

Abstract

Two methods are described to estimate the boundary of a 2-D penetrable object and the average sound speed in the object. One method is for circular objects centered in the coordinate system of the scattering observation. This method uses an orthogonal function expansion for the scattering. The other method is for noncircular, essentially convex objects. This method uses cross correlation to obtain time differences that determine a family of parabolas whose envelope is the boundary of the object. A curve-fitting method and a phase-based method are described to estimate and correct the offset of an uncentered radial or elliptical object. A method based on the extinction theorem is described to estimate absorption in the object. The methods are applied to calculated scattering from a circular object with an offset and to measured scattering from an offset noncircular object. The results show that the estimated boundaries, sound speeds, and absorption slopes agree very well with independently measured or true values when the assumptions of the methods are reasonably satisfied.

I. INTRODUCTION

Ultrasonic imaging for medical applications has the advantages of being non-ionizing at intensities used for diagnostic imaging, showing motion in real time, and being inexpensive in comparison to other imaging modalities such as magnetic resonance and X-ray tomographic imaging. Consequently, ultrasonic b-scan imaging is commonly used to detect and diagnose disease. Although b-scan imaging based on transmitting pulses and receiving echoes with focused beams is widespread, ultrasonic imaging based on inverse scattering is relatively undeveloped in medical applications. This lack of development can be attributed to several reasons. The most important of these is that a practical solution of multidimensional inverse problem does not yet exist and current solutions now require a priori knowledge or iteration or both. Other important reasons are the complex apparatus required to acquire data and the substantial computation needed to reconstruct images.

A common method for simplification of 3-d inverse scattering computations is to reduce the problem to a series of 2-D reconstructions. In this approach, the object is divided into cross sections that are isolated by transmitted fields and receiver sensitivity patterns that are focused in elevation. These scattering measurements form data sets that can be used in 2-D inversions of the cross sections. These 2-D problems are not as computationally intensive as the 3-D problem, but retain other difficulties that are inherent to inverse scattering computations of any dimension.

Because iterative reconstruction of large-scale and high-contrast objects must begin with a reasonably accurate initial estimate, at least a global knowledge such as object shape, position, and average sound speed, must be available before reconstruction of the object or region being imaged. A prior study [1] that used constraints to facilitate image reconstruction was aimed at finding the boundary of an obstacle. In that study, a linear sampling method [2] was used to obtain an a priori estimate of the obstacle boundary but the internal characteristics of the obstacle were not found. The purpose of this study is to examine ways in which the global characteristics of scattering objects in addition to the object boundary can be estimated. In the examination, emphasis is on scattering objects that are either radial or essentially radial because radial objects are a good model for applications such as imaging the breast. Both calculated scattering and measured scattering are used for the estimation of the boundary, average sound speed, and the offset of the object from a center of a coordinate system.

The choice of parameters in this study is motivated by an available ring transducer system [3] and imaging experiments [4], [5] conducted with the system using an eigenfunction method of inverse scattering [6]. In this system, the transducer is a 150-mm diameter ring of 2048 equally spaced transducer elements, each with an elevation of 25 mm. The elements are focused in elevation with an acoustic lens that creates a beam in the imaging plane with an elevation f-number of 5 in the center of the ring and which has a nominal center frequency of 2.5 MHz with an average −6-dB bandwidth of 70% around the center frequency. In the imaging experiments, 6-mm diameter objects with weak contrast and a 48-mm-diameter tissue-mimicking phantom were used. Although the weak contrast objects could be reconstructed using a Born approximation with a homogeneous background, the 48-mm phantom object was relatively large (about 100 wave-lengths at the center frequency) and had relatively high contrast (about 10%), so an initial estimate of the object characteristics and iteration were needed to reconstruct a satisfactory image.

The main idea of inverse scattering methods that use an estimate of the scattering object properties to obtain a detailed quantitative image from the roughly estimated scattering medium or, in other words, to reconstruct the scattering potential in a known large-scale and high-contrast inhomogeneous background is to linearize the nonlinear solution in a small range around the true solution. A requirement, however, is that the major features of the scattering medium be known to obtain a good approximation. In view of this, estimation of scattering object parameters such as shape and average sound speed is a fundamentally important first step in the reconstruction process.

Level-set methods [7] were originally used to find the boundary of impenetrable objects (i.e., obstacles) but these methods also have been adapted to acoustically penetrable objects [8]. These approaches characterize the interface between the background and the scattering medium by using the characteristics of the far-field scattering patterns from a large number of incident-wave directions. In this paper, a similar but less-complicated and more direct method is used to estimate the boundary of a penetrable object as well as other object characteristics such as average sound speed and attenuation.

The scattering object may also contain large-scale and high-contrast internal features that need to be included in the background to ensure convergence of inverse scattering calculations. Estimation of these features is not addressed by the methods discussed here, and is also not addressed by the level-set methods cited above. Consequently, development of additional techniques may be necessary to secure good starting points for iterative calculations. However, the influence of internal features in medically relevant applications such as cross-sectional reconstructions of the pendant breast is likely to be smaller and more localized than the influence of the boundary and average medium properties of the entire cross section. Thus, initial emphasis on recovering global characteristics is appropriate.

The use of a priori knowledge to begin the reconstruction process for general penetrable objects merits special comment. The scattering object boundary, the average sound speed in the object, and the offset of the object can all be derived from the same waveforms already needed for the reconstruction process. Thus, no additional acquisition of data is necessary. An initial average sound speed for the model that is close to the average sound speed in the true object lessens the need for iteration and, thereby, reduces the amount of computation as well as any risk that the iteration falls into a local minimum. For radial objects, the amount of computation is further reduced and the complexity of computation is reduced as well by centering the apparent origin of the measured data. The lessened amount of computation is another reason that estimation of the scattering object parameters is important.

The remainder of this paper is comprised of the following sections. Section II presents the basic relations and methods used for estimation of scattering object characteristics. Results for calculated and measured scattering are presented in Section III. Section IV discusses key assumptions and other steps used in the analysis and computations, and also discusses the results. Conclusions are presented in section V.

II. BASIC RELATIONS AND METHODS

A. Scattering Calculations and Measurements

For estimation of scattering object parameters from calculated scattering, a cylindrical object and an ellipsoidal object were examined. The calculation for a cylindrical object used an exact solution for a nonradial object consisting of multiple arbitrary-diameter nonconcentric cylinders with parallel axes [9]. The calculation for an elliptical object used a k-space method in which spatial derivatives are implemented in the spatial-frequency domain, the time history is obtained in steps using finite differences, and a perfectly matched boundary is included to avoid reflections from the edge of the spatial domain of the computation [10].

The analysis for a cylindrical object begins with consideration of a fluid cylinder that has radius acyl, sound speed ccyl, and density ρcyl and is positioned with its axis along the z-axis in a Cartesian coordinate system in which the coordinates in the xy plane are denoted (x, y) and the corresponding coordinates in a cylindrical coordinate system are denoted (r, ϕ). The cylinder is assumed to be illuminated by a monochromatic plane wave that has time dependence ejft and is traveling in a direction defined by a unit vector u(ϕi) in the xy plane. By matching the conditions of continuous pressure and normal component of particle velocity at the boundary of the cylinder, the pressure scattered from this cylinder in a lossless fluid background in which the sound speed is c and density is ρ can be expressed using an orthogonal function expansion in cylindrical coordinates as [11], [12]

ps(r,ϕ,ϕi)=m=+jmDm(f,acyl,ccyl)Hm(1)(kr)ejm(ϕϕi),    r>acyl.
(1)

In this expression, r is the radial distance in the xy plane, ϕ is the azimuthal angle of the observation point, ϕi is the azimuthal angle of the incident wave, k is the wavenumber of the incident wave, kcyl is the wavenumber in the cylinder, Dm are coefficients depending on the conditions at the boundary of the cylinder, and Hm(1) are cylindrical Hankel functions of the first kind. The coefficients Dm in the expansion are given by

Dm(f,acyl,ccyl)=Jm(kcylacyl)Jm(kacyl)γJm(kcylacyl)Jm(kacyl)γJm(kcylacyl)Hm(1)(kacyl)Jm(kcylacyl)Hm(1)(kacyl),
(2)

in which γ = (cρ)/(ccylρcyl) and ′ denotes differentiation with respect to the argument. The wavenumber k in the background and the wavenumber kcyl in the cylinder are given by

k=2πf/c,
(3)

and

kcyl=2πf/ccyl,
(4)

respectively.

The expression for Dm in (2) is also valid for a lossy cylinder with absorption coefficient αcyl, in which case the wavenumber kcyl is given by

kcyl=2πf/ccyl+jαcyl.
(5)

The incident pulse waveform g(t) in the calculations was a sinusoid with Gaussian envelope. The center frequency f0 was 2.5 MHz and the −6-dB bandwidth was 1.5 MHz. An explicit expression for the pulse is

g(t)=Aet2/2σt2sin(2πf0t),
(6)

in which the pulse-width parameter σt and −6-dB bandwidth B are related by

σt=2ln2πB
(7)

and A is a constant. The sound speed c and the density ρ in the background were 1.509 mm/µs and 0.997 g/cm3, respectively. The background was lossless.

For calculations of scattering from a cylindrical object, the diameter of the cylinder was 48.97491 mm; the center of the cylinder was offset from the origin of the coordinate system of the calculation. The offset angle was 183.8529° and the offset distance was 0.2615 mm. For calculations of scattering from an elliptical object, the semimajor axis and the semi-minor axis were 25.1296 mm and 23.8758 mm, respectively. The semi-major axis was inclined 1.5366° above the horizontal axis and the center of the ellipse was offset by 0.2695 mm at an angle of 180.3189°. Each object enclosed two 10-mm (nominal) diameter cylinders and three 0.1000-mm diameter cylinders. The larger enclosed cylinders had respective diameters of 10.5418 and 10.5082 mm and centers at (9.0572, 2.0106) and (−9.6678, −1.6991) mm. The smaller enclosed cylinders had had centers at (6.5354, 12.9925), (−11.4739, 9.3690), and (−7.4163, −12.6167) mm. The large and small enclosed cylinders had respective sound speeds of 1.576 and 2.760 mm/µs, densities of 0.940 and 1.580 g/cm3, and absorption-coefficient slopes of 0.341 and 2.797 dB/(cm·MHz). In the enclosing cylindrical object and in the enclosing ellipsoidal object, the sound speed was 1.574 mm/µs, the density was 1.040 g/cm3, and the absorption-coefficient slope was 0.458 dB/(cm·MHz). For both the cylindrical and ellipsoidal objects, the number of incident wave directions, i.e., views, was 1024, the number of observation positions was also 1024, the temporal sampling rate was 10 MHz, and the time interval was 0.10 µs. The angles of the views and the observations were the same and evenly distributed around the unit circle. All the parameters for the calculations were chosen so that the calculations of scattering mimic the measurements of scattering.

For estimation of scattering object parameters from measured scattering, the incident and scattered waves were measured using the previously noted unique ring transducer system developed for scattering and imaging studies. The data acquisition for both the total field and incident field spanned several hours because the measurement system synthesized the transmit aperture that implemented the spatially limited incident plane wave and also synthesized the receive aperture that recorded the total pressure around the circumference of the measurement ring [3], [5]. During the measurement period, a temperature control system limited the temperature fluctuations in the measurement chamber to less than ±0.2°C. As detailed in [5], the incident wave was subtracted from the total wave after the waveforms were compensated for time shifts caused by temperature fluctuation by using a reference beam that propagated across a secant of the measurement ring without intersecting the scattering object. Figures showing the total field (before subtraction) and scattered field (after subtraction) are in [5, Fig. 5]. A reference waveform was calculated by backpropagating the measurements of the incident wave to the center of the ring and averaging the resulting waveforms for each view laterally across a 50-mm interval around the center of the ring. The resulting scattered pressure was used for estimation of the scattering object parameters.

Fig. 5
Boundary estimations. In the left panel, the solid red line and the solid blue line denote the boundaries derived from calculated scattering by an offset cylindrical object and by an offset elliptical object, respectively. The corresponding color dashed ...

The scattering object used for measurements was a cylindrical tissue -mimicking phantom with a nominal 48-mm diameter, an average sound speed of 1.574 mm/µs, an average density of 1.040 g/cm3, and an average attenuation coefficient with a slope of 0.458 dB/(cm·MHz) at 29°C. (Because attenuation is the sum of absorption loss and scattering loss and the scattering loss in the phantom is known to be small compared with the absorption loss, the attenuation coefficient slope was also considered to be the absorption slope of the absorption coefficient.)

B. Object Parameter Estimation

1) Boundary and Sound Speed Estimation for Cylinders

A monochromatic incident plane wave e jku(ϕir traveling in a direction ϕi and illuminating a cylinder of radius acyl with associated wavenumber kcyl results in a scattered wave at the observation ring given by

φscattheory(rring,ϕϕi,f,acyl,ccyl)=m=+Dm(f,acyl,ccyl)jmejm(ϕϕi)Hm(1)(krring),
(8)

where the Dm(f, acyl, ccyl) are given by (2). The product Dm(f,acyl,ccyl)jmHm(1)(krring) in the above equation can be interpreted as the mth coefficient

φ^scattheory(rring,m,f,acyl,ccyl),
(9)

in the Fourier series expansion

φscattheory(rring,ϕϕi,f,acyl,ccyl)=m=+φ^scattheory(rring,m,f,acyl,ccyl)ejm(ϕϕi).
(10)

This expansion can also be formed for the measured data

φscatmeas(rring,ϕϕi,f)=m=+φ^scatmeas(rring,m,f)ejm(ϕϕi)
(11)

by taking a discrete Fourier transform of an appropriately normalized view in the angular direction.

A natural way to assign the model parameters acyl and ccyl is to fit the theoretical coefficients

φ^scattheory(rring,m,f,acyl,ccyl)=Dm(f,acyl,ccyl)jmHm(1)(krring)
(12)

to the measured coefficients

φ^scatmeas(rring,m,f).
(13)

In other words, adjust acyl and ccyl to get the best fit for

Dm(f,acyl,ccyl)~φ^scatmeas(rring,m,f)jmHm(1)(krring),
(14)

in which the left side is determined by the theoretical model and the right side is obtained from the measurements. However, the number of measurements available to fit is very large, i.e.,

(Num of Meas)=(Num of Fourier Coef)×(Num of Freq)×(Num of Views),
(15)

where Num, Meas, Coef, and Freq are abbreviations for number, measurements, coefficients, and frequencies, respectively. This suggests that fitting be performed in a least-squares sense and need only involve a subset of values in the above count. The details of the fitting process should be based on an examination of the behavior of dependency of Dm(f, acyl, ccyl) with respect to each of its arguments.

Cross correlation can also be used to estimate the position of a scattering object boundary and the average speed of sound in the scattering object. Cross correlation can be used in addition to find the offset of the scattering object from the coordinate system origin. In the cross-correlation approach, a window is applied to isolate the waveform from the region of interest and the waveform is cross correlated with a reference waveform. The position of the cross-correlation peak is then taken to be the arrival time of the scattered wave.

The cross-correlation approach is based on relations developed from consideration of a plane wave propagating at an azimuthal angle ϕ and impinging on an object that scatters the wave. The relations are obtained by assuming that the wave backscattered at angle ϕ − 180° is observed at a point E on a ring of radius R, that the center of the observation ring is the origin of the spatial coordinates, that the time origin is chosen so a reflection from the center of the observation ring arrives at time t = 0, and that t180° is the arrival time of the first reflection from the boundary of the object. This choice implies that the reference pulse is transmitted at time t = −2R/c. Geometry considered to develop the relations is shown in Fig. 1.

Fig. 1
Geometric configuration for boundary estimation. A plane wave denoted by parallel dotted lines and an accompanying arrow is incident on a scattering object denoted as a shaded region with an irregular boundary. Backscattered reflections from the scattering ...

If the first reflection occurs at a point on a ray between the observation point E = Re(ϕ) and the center of the ring, then the coordinate of that point will be r(ϕ)e(ϕ), where r(ϕ) is obtained from t180° by using r(ϕ) = −ct180°/2 and e(ϕ) is a unit vector in the direction of the observation point E. However, reflections from other locations can also arrive at time t180° so the location of the first reflection is not uniquely determined. For example, if the reflection point is offset a distance s from the ray at angle ϕ, then the new path length rpar(s) from the offset reflection point to the same observation point E will be longer than the original path length rpar(0) = Rr(ϕ) for a reflector at r(ϕ)e(ϕ) and the distance that the incident wave propagates to the offset reflection point must be shorter to maintain the same overall transmission and reflection distance. The total propagation distance in the case of an offset s will be the same as in the case of the intersection along the ray with no offset if the difference between the offset-reflection path length rpar(s) and the reduction rred(s) in the original path length is equal to the difference between the radius of the observation ring and the distance to the boundary in the no-offset case, i.e., rpar(s) − rred(ϕ) = Rr(ϕ). The locus of points [ell](s) satisfying this condition forms a parabola that can be expressed as

(s,ϕ)=se(ϕ)+{s24[Rr(ϕ)]+r(ϕ)}e(ϕ),
(16)

in which e[perpendicular](ϕ) is a unit vector orthogonal to e(ϕ) in the counter-clockwise direction and the angular dependence of e, e[perpendicular], and r on the azimuthal angle ϕ is explicitly indicated. (This motivates the use of the subscript, used to designate the earliest reflection distance.) As the azimuthal angle of the incident wave ranges from −180° to +180°, the locus of reflection points corresponding to the earliest reflection distance r(ϕ) forms a one-parameter family of parabolas. Choosing the earliest reflection for each incident wave angle insures that the parabola [ell](s, ϕ) is tangent to the boundary of the scattering object at the point of intersection.

Let s(ϕ) denote the intercept parameter for the parabola [ell](s, ϕ). Then, using (16), a parametric expression g(ϕ) for the boundary of the scattering object is

g(ϕ)(s(ϕ),ϕ)=s(ϕ)e(ϕ)+{s2(ϕ)4[Rr(ϕ)]+r(ϕ)}e(ϕ).
(17)

In this expression, the function s(ϕ) is found by imposing the constraint that g(ϕ) be tangent to the parabola [ell](s, ϕ). This is equivalent to requiring that the tangent of g(ϕ) be orthogonal to the normal of [ell](s, ϕ) at their intersection.

At the point s(ϕ) where the parabola and boundary intersect, the angular derivative relations (d/dϕ)e(ϕ) = e[perpendicular] (ϕ) and (d/dϕ)e[perpendicular] (ϕ) = −e(ϕ) can be used to express the tangent Tϕ to the parabola [ell](s, ϕ) in the form

Tϕ=s(s(ϕ),ϕ)=e(ϕ)+s(ϕ)2[Rr(ϕ)]e(ϕ).
(18)

The corresponding form for the normal Nϕ to the parabola [ell](s, ϕ) at the point s(ϕ) is

Nϕ=e(ϕ)+s(ϕ)2[Rr(ϕ)]e(ϕ).
(19)

With this form for the normal Nϕ, the condition that g(ϕ) be tangent to the parabola [ell](s, ϕ) can be expressed using an inner product as the relation

0=Nϕ,ddϕg(ϕ).
(20)

Because the partial derivative of [ell](s, ϕ) with respect to s is orthogonal to the vector Nϕ at the point of intersection, this equation reduces to

0=Nϕ,ϕ(s(ϕ),ϕ).
(21)

The angular derivative relations can then be used again to obtain

ϕ(s(ϕ),ϕ)=s(ϕ)e(ϕ)+{1+s2(ϕ)4[Rr(ϕ)]2}r(ϕ)e(ϕ)+{r(ϕ)+s2(ϕ)4[Rr(ϕ)]}e(ϕ).
(22)

Substituting this result in (21) and rearranging terms yields

0=Nϕ,ϕ(s(ϕ),ϕ)=r(ϕ)+{1+r(ϕ)2[Rr(ϕ)]}s(ϕ)+{r(ϕ)4[Rr(ϕ)]}s2(ϕ)+{18[Rr(ϕ)]2}s3(ϕ).
(23)

For nearly radial objects, the quadratic and cubic terms in s are negligible. Then, solving the resulting linear equation for s(ϕ) gives

s(ϕ)=r(ϕ)/{1+r(ϕ)2[Rr(ϕ)]}.
(24)

Eqs. (17) and (24) together determine the boundary of the scattering object based on the earliest reflection times in the backscatter direction.

If the arrival time of a reflection from the center of the observation ring cannot be determined, the distance r(ϕ) cannot be obtained by using r(ϕ) = −ct180°/2. However, an alternative method based on time differences can be used when the scattering object is nearly radial and first-order scattering dominates. Then, the arrival time of the earliest backscattered wave can be considered to be the arrival time of the echo of the nearest boundary of the object and the arrival time of the latest backscattered wave can be considered to be the echo from the farthest boundary of the object. The time difference between the positions of these two arrival times is the time of propagation through twice the dimension of the object in the direction of transmission. This difference Δt180°(ϕ) in the backscatter direction at angle ϕ can be expressed

Δt180°(ϕ)=2dobj(ϕ)/c¯obj,
(25)

in which cobj is the average speed of sound inside the object along the length of the dimension being traversed and dobj(ϕ) is the dimension of the object in the direction defined by the angle ϕ. The positions of the peaks are found, as described earlier, by cross correlation of the backscattered waveform with a reference waveform.

Also, for nearly circular scattering objects and first-order scattering, the time difference between the positions of first and last peaks in the forward-scattered wave is the difference between the time of propagation through the dimension of the object in the forward scattering direction and the time of propagation through the same path in the background medium. This time difference Δt in the forward scattering direction can be expressed

Δt0°(ϕ)=dobj(ϕ)(1/c1/c¯obj),
(26)

where c is, again, the sound speed in the background medium. The positions of these peaks are conveniently found by cross correlation of the forward-scattered waveform with a reference waveform.

The backscatter time difference Δt180° given by (25) and the forward scatter time difference Δt given by (26) can be used to calculate the average sound speed cobj and the dimension dobj of the scattering object in every transmit direction. The estimate of the average sound speed in the object is then found to be

c¯obj=[1+2Δt0°(ϕ)/Δt180°(ϕ)]c,
(27)

and the estimate of the dimension of the object is found to be

dobj=c[Δt180°(ϕ)/2+Δt0°(ϕ)].
(28)

These estimates of the diameter and sound speed of objects only use the forward scattering and backscattering measurements in a single view but can be computed for each view. Because the incident wave angle ϕ and the incident wave angle ϕ − 180° are along the same line but in opposite directions, the values of Δt180°(ϕ) and Δt180°(ϕ − 180°) represent the same time difference and the values dobj(ϕ) and dobj(ϕ − 180°) represent the same dimension. Thus, the average values [cobj(ϕ) + cobj(ϕ − 180°)]/2 and [dobj(ϕ) + dobj(ϕ − 180°)]/2 can be used as the average sound speed in the object and the dimension of the object along the line defined by ϕ. If the object is radial, even better estimates of cobj and dobj can be found by averaging the cobj and dobj values obtained from (27) and (28) for the views with ϕ ranging over a 180° interval.

As the shape of the scattering object deviates from a perfect cylinder, (28) loses accuracy as an estimate for the diameter of the object along the e(ϕ)-axis. However, the object boundary can still be completely determined as the envelope of a family of parabolas as described earlier. This computation employs the distance function r(ϕ) = (c/2)t180°(ϕ), where t180°(ϕ) denotes the arrival time of the earliest backscattered echo in the e(ϕ)direction. The times Δt180°(ϕ) and Δt(ϕ) that appear in (25), (26), (27), and (28) are differences and, therefore, independent of the time origin. However, t180°(ϕ) must be measured on a time scale in which a reflection from the center of the ring arrives at time t = 0. This time origin can be established from measurements of the unobstructed pulse wave because back-propagation of these measurements allows identification of the time that the pulse arrives at the center of the ring. The distance function r(ϕ) is used in (24) to determine s(ϕ). Then, s(ϕ) is used in (17) to determine the boundary γ(ϕ). If the object is circular and centered on the center of the ring then the equation for the boundary is just γ(ϕ) = r(ϕ)e(ϕ), where r(ϕ) is constant.

The positions of correlation peaks were found in this study by finding the maximum of the spline-interpolated cross -correlation function in the neighborhood of the maximum computed value. The interpolation of the cross-correlation function improves the resolution of the position. However, a small difference between c and cobj can make identification of correlation peaks unreliable when only one waveform is used. Averaging multiple waveforms can improve peak identification. An average of five waveforms (two on each side of the observation position and the waveform at the observation position) was used in this study to estimate the reflection and transmission times in the calculation of cobj and dobj.

The boundary γ (ϕ) of the scattering object can be represented by (17) and the term s(ϕ) appearing in (17) can be approximated by (24) for objects of any shape. However, approximations leading to (24) become less accurate for irregular or eccentric objects. Although the argument ϕ in the boundary curve of an irregular object no longer represents azimuthal angle, the boundary can, nevertheless, be converted to polar form by interpolating the discretely sampled boundary vectors {γ (ϕ1), …, γ (ϕN)} to cylindrical coordinates.

2) Absorption Estimation for Cylinders

The total power loss Ptot from a monochromatic plane wave incident on a scattering object can be decomposed into a scattering component Pscat and an absorption component Pabs. This can be expressed

Ptot=Pscat+Pabs.
(29)

An expression for the total power loss is given in terms of the normalized far-field pattern of scattering A(θ)[6] by the extinction theorem2 [14], [15]

12ωρ0Im[A(0)]=k02ωρ0Im[x2+y2<Bdryη(x,y)|p(x,y)|2dxdy]+1412ωρ0[12π02π|A(θ)|2dθ]
(30)

in which the left side is the total power loss, the first term on the right side is the power loss caused by absorption, the second term on the right side is the power lost because of scattering, k0 is the wavenumber in the background medium, η is the complex-valued index of refraction in the scattering object, and p is the total pressure. Values of Pabs can then be computed as differences: Pabs = PtotPscat for either calculated or measured data. Because estimates of sound speed and the boundary are known, the absorption slope β can be varied to fit Pabs determined as a function of temporal frequency. Because the scattering loss in the scattering object used for measurements is, as noted earlier, known to be small compared with the absorption loss and because the choice of parameters for the calculation of scattering was made to mimic scattering by the object used for measurements, the absorption coefficient slope in this study can also be considered to be the slope of the attenuation coefficient.

C. Scattering Object Offset Estimation and Correction

Estimation and correction of the object offset in measured (or calculated) scattering data can simplify calculations that accomplish the reconstruction of the object, particularly when the object is circular (see, e.g., [5]). If the scattering object is circular or elliptical or has symmetry, then the offset of the center of the object from the origin of the observation ring can be estimated and corrected by using a so-called geometric method or a so-called phase method.

1) Geometric Method

A circle and an ellipse were used to fit the estimated boundary of the object in this study because the measured scattering was from a tissue-mimicking phantom that is nominally a cylinder but has a cross section that is nearly an ellipse with a small eccentricity. Least-mean-square-error methods were used to determine the center and radius of the fitted circle and also to determine the foci, major axis, and minor axis of the fitted ellipse. The offset of the center of the fitted circle or ellipse was considered to be the offset of the object.

2) Phase Method

A representative multiple scattering path produced by a monochromatic incident plane wave is illustrated in Fig. 2 for a scattering object and the same object shifted by the small amount δ. In this figure, each arrow is associated with a transmission factor that depends on the length of the arrow and the arrow endpoints are associated with reflection coefficients that are determined by the medium properties at that location. The total scattered wave is a (continuous) sum of such contributions. The figure also shows unit vectors that are used to describe the offset.

Fig. 2
Multiple scattering in an unshifted and shifted object. The boundary of the unshifted object is shown as a closed irregular contour. The boundary of the object shifted by δ is shown by a corresponding dotted contour. An incident wave traveling ...

For both the unshifted and shifted objects, the internal transmission and reflection factors are the same so the difference between the scattering produced by the two paths is entirely due to the changes in length of the initial and final arrows. This effect is accurately approximated by the factor

ejk(uius)·δ
(31)

in which ui is a unit vector in the direction of the incident wave and us is a unit vector in the direction of the observation point. Because this factor is the same for every displaced path, translation of the scattering volume causes the total scattering measurement at the transducer element at angle us to be multiplied by the phase factor e jk(uiusδ, i.e.,

φscat(δ)(ui,us)=φscat(0)(ui,us)ejk(uius)·δ.
(32)

Thus, offset-corrected scattering data can be obtained using

φscat(0)(ui,us)=φscat(δ)(ui,us)ejk(uius)·δ.
(33)

The pair of vectors ui, us is conveniently parameterized for analysis by the angle θ between them.

The problem is to determine the offset δ from the scattering measurements. An object of general shape does not, however, have a preferred center so no way exists to make this assignment. Nevertheless, if the object has symmetries, then the problem becomes meaningful. For example, suppose that the scattering object is symmetric with respect to reflection about the axis of the incident beam so that

φscat(0)(ui,us(θ))=φscat(0)(ui,us(θ)).
(34)

In this case, the product (of measured values)

φscat(δ)(ui,us(θ))φscat(δ)(ui,us(θ))¯  =φscat(0)(ui,us(θ))φscat(0)(ui,us(θ))¯×ejk[uius(θ)]·δejk[uius(θ)]·δ  =|φscat(0)(ui,us(θ))|2ejk[us(θ)us(θ)]·δ
(35)

has phase angle

Angle[φscat(δ)(ui,us(θ))φscat(δ)(ui,us(θ))¯]=k[us(θ)us(θ)]·δ.
(36)

Reference to Fig. 2 indicates that this angle can be expressed

Angle[φscat(δ)(ui,us(θ))φscat(δ)(ui,us(θ))¯]=2kδsinθ
(37)

in which δ[perpendicular] is the component of δ orthogonal to the direction ui of the incident wave.

The result given in (37) shows that uncomplicated calculations involving the scattering measurements for a single view and at a single frequency permit evaluation of the component of the offset vector δ along an axis that is orthogonal to the incident beam. Thus, scattering at a single frequency f in two orthogonal views is sufficient to completely determine the offset. The redundancy provided by 1024 views and multiple frequencies can be employed to reduce the statistical error from noise in this determination. This calculation should be carried out in advance of the fitting process used to determine the radius and sound speed of the model cylinder and the factor ejk(uiusδ should be removed from the scattering measurements before modeling.

Experience with measured data indicates, however, that phase values in the presence of noise are inaccurate for small angles. This inaccuracy can be overcome by computing phase differences for larger angles. However, to account for phase wraps, the phase at larger angles must be computed by summing phase differences in small increments of θ starting with a known value. The basic relation is

Phase(θ0+NΔθ)=Phase(θ0)+[Phase(θ0+Δ)Phase(θ0)]+[Phase(θ0+2Δ)Phase(θ0+Δ)]++[Phase(θ0+NΔ)Phase(θ0+(N1)Δ)],
(38)

or, equivalently,

Phase(θ0+Δθ)=Phase(θ0)+AccumDiff().
(39)

Using (38), the accumulated phase differences, Accumdiff, can be written

AccumDiff()=2kδsin(θ0+Δθ)Phase(θ0),=1,2,,N.
(40)

This expression has the form of a linear fit with the two unknown coefficients δ[perpendicular] and Phase(θ0) and can be written in the matrix form

(AccumDiff(0)AccumDiff(1)AccumDiff(N))=(ksin(θ0)1ksin(θ0+Δθ)1ksin(θ0+NΔθ)1)(2δPhase(θ0)).
(41)

The minimum-mean-square-error solution for δ[perpendicular] and Phase(θ0) can, therefore, be expressed

(2δPhase(θ0))=(QtQ)1Qt(AccumDiff(0)AccumDiff(1)AccumDiff(N)),
(42)

in which the matrix Q is given by

Q=(ksin(θ0)1ksin(θ0+Δθ)1ksin(θ0+NΔθ)1).
(43)

An alternative method also based on the symmetry of the scattered wave around the forward (0°) direction uses a variational approach to estimate the offset of a radial object. In this method, the phase differences on the left side of (37) are fitted to the analytic expression on the right side of the equation by a search through a range of δ[perpendicular] values. The component of δ orthogonal to δ[perpendicular] is found by a similar fit to the scattering from an orthogonal view.

The steps in the variation-of-parameters method to find δ[perpendicular] are:

  1. For a given δ[perpendicular], scattering angle, and temporal frequency, compute an exponential phase factor corresponding to the analytic expression on the right side of (37) and multiply the product of the fields on the left side of (37) by the conjugate of this factor to find the residual phase of the offset.
  2. Evaluate the residual phase for a range of scattering angles and frequencies and calculate the mean-square value of these residual phases.
  3. Repeat step 1 and step 2 above for different values of δ[perpendicular] to find the minimum-mean-square value of the residual phase.

III. Results

A. Scattering Calculations and Measurements

Representative waveforms and cross correlations for calculated scattering and for measured scattering are shown in Fig. 3. The reference waveform was obtained by first backpropagating to the center of the ring measurements of the unobstructed incident wave and then averaging the waveforms across the spatially uniform extent of the wave as described in [5]. The temporal positions of the cross-correlation peaks in the backscatter direction define the arrival time of the near-wall echo and the far-wall echo from the scattering object.

Fig. 3
Representative waveforms and cross correlations for arrival time estimation from calculated and measured scattering. The dashed lines are for calculated scattering and the solid lines are for measured scattering.

B. Estimations of Scattering Object Parameters

The radius r and sound speed c estimated using the frequency and index dependence of Dm to obtain a minimum-mean-square-error fit to (14) are shown in Table I for calculated scattering from a centered circular object, an offset circular object, and for measured scattering from a nearly circular offset object both before and after offset correction. The fit used values of m and f in the ranges of 75 to 125 and 2.50 to 3.25 MHz, respectively, because calculations of the Dm dependency on index and frequency indicated that the Dm have regular peaks as a function of index and vary smoothly as a function of frequency in these ranges for the object from which scattering was measured. The results show that the Dm method works well for a centered circular object because the error between the estimated and true radius and the estimated and true sound speed is at most a fraction of a percent. However, the results also show parameter estimations using the Dm method for an offset circular object are not as accurate because the error in radius and in sound speed increases to more than 10 times the corresponding errors for a centered object with a relatively small offset of 0.2615 mm. Thus, the Dm method only works for objects that are almost perfectly centered, circular, and relatively homogeneous.

TABLE I
True and Estimated Values of Radius and Sound Speed Found Using the DM Method for Calculated Scattering.

The average sound speed and dimension estimated using (27) and (28) for calculated scattering from an offset circular object, from an offset elliptical object, and for measured scattering from an offset noncircular object are plotted in Fig. 4 as a function of view, i.e., angle around the scattering object. The average sound speed averaged over views is 1.5751 mm/µs and 1.5764 mm/µs for the off-set circular object and the offset elliptical object, respectively. These values compare favorably to the true value of 1.5740 mm/µs and indicate that a small offset of the scattering object from the center of the coordinate system has little effect on the calculation of the backscatter time difference given by (25) and the forward scatter time difference given by (26).

Fig. 4
Values of object parameters, i.e., dimension d and average sound speed c, estimated using the cross-correlation method and the parabola method in combination with calculated scattering from an offset circular object, from an offset elliptical ...

The boundary obtained from using (17) and (24) for calculated scattering from an offset circular object and for an offset elliptical object and for measured scattering from a noncircular offset object are plotted in Fig. 5. The root-mean-square error (rmse) of boundaries between the estimated values and the true values for circle and ellipse are 0.0039 mm and 0.0048 mm, respectively. The small errors show that the cross-correlation and parabola methods used to estimate the boundary in these cases are accurate. Application of the methods to measured data yields a radius of 24.4875 mm and an offset distance and angle of 0.2615 mm and 183.8529°, respectively, for the fitted circle and foci of f1 = (7.5662, 0.2087) mm, f2 = (−8.1053, −0.2117) mm and a semi-major axis of 25.1296 mm and a semi-minor axis of 23.8758 mm for the fitted ellipse. The rmse of the fitted circle and the fitted ellipse is 0.4881 mm and 0.2069 mm, respectively. The smaller error for the fitted ellipse shows that the scattering object used for measurements is more like an ellipse with a small eccentricity than a than a circle. The average sound speed estimated from the measured data are 1.5724 mm/µs. Boundary estimations made from calculations of scattering from centered objects were indistinguishable from the true boundaries.

Averages of absorption coefficients estimated from measured scattering data by using the extinction theorem are plotted in Fig. 6 for three sets of views along with linear approximations fitted to the averages as a function of frequency. The slope fitted to the average of all the views is 0.459 dB/(cm·MHz), a value that closely agrees with the true value of 0.458 dB/(cm·MHz). The slope fitted to average of views with consecutive illumination of the internal cylinders is less than the slope fitted to the average of views with separate illumination of the internal cylinders because the internal cylinders have less absorption than the cylinder that encloses them.

Fig. 6
Absorption coefficient. The solid lines are averages of absorption coefficients calculated for three different sets of measured-scattering views by using the 2-D version of the extinction theorem. The dashed lines are linear least-mean-square-error fits ...

C. Object Offset Estimation

The offsets obtained from the fit of a circle and the fit of an ellipse to calculated scattering and measured scattering are given in the foregoing paragraph. Offsets found by applying the phase method to calculated and measured scattering are shown in Table II. The results show the phase method performs well with calculated scattering in the absence of noise. The results also show that offsets from the phase method using measured scattering are corrupted by noise and that values found using parameter variations are relatively insensitive to noise but still unacceptable.

TABLE II
True and Estimated Values of Offset in a Direction Orthogonal to the Transmit View of an Offset Circular Scattering Object.

D. Offset Correction and Image Reconstruction

Images reconstructed from calculated scattering data and from measured scattering data after offset correction are shown in Fig. 7. True values of sound speed and attenuation slope and the corresponding initial estimates as well as the final values in the reconstructions are listed for comparison in Table III. Values in the final estimates column are an average in a (representative) 9.8-mm diameter circular region centered at (8.0, −10.0) mm in each image. The offset correction of calculated scattering data used (33) with the offsets in Table II and the offset correction of measured scattering data used a fitted circle. Image reconstruction employed the inverse scattering method described in [5] for high-contrast large-scale objects by using the foregoing estimations of the average sound speed, radius, and absorption as a priori knowledge to calculate scattering by the inhomogeneous background in the first reconstruction, i.e., the zeroth iteration. Because the initial estimates obtained using the methods described were accurate and the scattering object variations about the mean values were small except for three filaments that acted like point scatters, the reconstruction, which lacks speckle even though the phantom contains small scatters that produce speckle in b-scan image, converged during the next reconstruction, i.e., the first iteration.

Fig. 7
Image reconstructions. The images in the top row are from offset-corrected calculated scattering. The images in the bottom row are from offset-corrected measured scattering. The scattering object is a 48-mm diameter tissue-mimicking phantom. The phantom ...
TABLE III
True or Measured and Estimated Sound Speed and Absorption Slope of the Scattering Objects.

IV. Discussion

The results of the Dm method show that the method can provide satisfactory estimates of sound speed and radius for a centered circular object but that the method may not provide satisfactory estimates of average sound speed and radius for a offset object. The reason for poor performance with an offset object is that the radius appearing in the argument of the Hankel function Hm(1) used to calculate Dm in (12) is relative to the center of the object rather than the center of the observation ring. This offset between the two origins must be corrected to obtain Dm that represent the scattering accurately. Therefore, compared with the Dm method, the cross-correlation method can provide more satisfactory results for boundary and sound speed estimations.

The attribution of arrival-time differences to offset and to noncircularity merits comment. If the scattering is from a cylindrical object perfectly centered in the observation ring, the time difference in the backscatter direction is constant. If the cylindrical object has a 48.9749-mm diameter and is offset 0.2615 mm, the maximum backscatter time difference is found to be less than 1 ns. However, the corresponding maximum time difference found from the measured scattering by a similar size object is about 3 µs. Because the offset of the object from which scattering was measured is known to be small (i.e., about a quarter of a millimeter), this large backscatter time difference is attributed to noncircularity of the scattering object.

The offset-estimation accuracy in the vertical direction is not as good as in the horizontal direction. This is because the method is based on the approximation that the time difference in backscatter direction can be used to represent the dimension of the object in every transmit direction through the center of ring. This approximation is valid only when the offset distance is very small, i.e., less than 1 mm, and the accuracy of the offset estimation decreases as the offset distance increases. Because the horizontal offset of the considered objects is bigger than the vertical offset, the estimation accuracy in the vertical dimension is poorer than in the horizontal dimension. Nevertheless, the error in the estimates of offset in the current study is acceptable.

The importance of correcting measured scattering data for the offset in the center of a nearly circular object from the center of the observation ring when the reconstruction process uses eigenfunctions of the scattering operator, as is the case for the reconstructions shown in Fig. 7, is noteworthy. As detailed in [5], the number of calculations to obtain fields used in the reconstruction is reduced by about an order of magnitude by symmetries and the availability of analytic expressions for scattering by the homogeneous model object used in the reconstruction process. Use of measured scattering that has not been corrected for offset with model scattering data that is from a centered object results in prominent crescent artifacts at the boundary of the reconstructed scattering object. This artifact is notably absent in the image reconstruction from offset-corrected calculated scattering. The residual artifacts on the boundary of the reconstructions using offset-corrected measured data result from the noncircularity of the scattering object when circularly symmetric model is used in the reconstruction process.

The images reconstructed from calculated scattering have very clear boundaries and internal features inside as shown in Fig. 7. As also shown in Fig. 7, the images reconstructed from the measured scattering are very good, although some artifacts are present. The artifacts are attributed mainly to three causes. The first is the inevitable crosstalk in the measurement system. The second is that the scattering object is not exactly circular. The third is the incomplete removal of the incident wave from the total wave.

The procedure by which images like those in Fig. 7 are reconstructed for large-scale high-contrast objects merits special comment. Results presented in [5] show that reconstructions of such objects do not converge to the true values of sound speed and absorption slope when the reconstruction process begins with an empty background. Thus, good initial estimates of object characteristics are important not only to reduce iteration but also to ensure converge of iterative reconstruction processes such as described in [5] for large-scale, high-contrast objects. The accuracies of the final sound speed and absorption slope estimates reported in Table III are influenced by a variety of factors, such as the number of iterations and the resolution of the reconstruction algorithm, that are unrelated to the accuracy of the initial estimates. However, the final estimates in Table III demonstrate that enough precision was secured from the initial estimates to ensure convergence of the iterative reconstructions. This is the criterion used to determine that the a priori estimates are satisfactory.

The method that finds the boundary of a scattering object using the envelope of a family of parabolas was investigated for objects with elliptical boundaries of varying size and eccentricity and also for objects with boundaries having radii of the form r([var phi]) = r0 + Δr cos (N [var phi]), where Δr and N are, respectively, the amplitude and spatial-frequency parameters that modulate the boundary. The reconstructions, as expected from theoretical considerations, were essentially exact for the elliptical boundaries. However, the reconstructions underestimated depressions in the modulated boundaries when the curvatures of the depressions exceed the curvatures of the parabolas. The validity of the approximation in (24) was also investigated for these examples and was found to be highly accurate for boundaries in geometries such as those anticipated for coronal sections of a pendant breast in medical applications of ultrasound.

V. Conclusion

Two methods have been developed to estimate the boundary of a scattering object and the average sound speed in the object. One method is based on a cylindrical function expansion of the scattering and is limited to circular objects that are centered in the observation coordinate system. However, curve-fitting and phase-based methods are described to estimate and correct observed scattering so the method can also be applied to offset circular objects. Another method uses cross correlation of received waveforms and represents the object boundary as the envelope of parabolas. This method is applicable to noncircular objects with an essentially convex boundary. Application of the methods to calculated scattering by cylindrical and elliptical objects with and without an offset, and also to measured scattering from an offset noncircular object illustrates the effectiveness of the methods in realistic situations. The offset estimation method that fits a circular or elliptical contour to a boundary is limited to objects that are nearly circular or elliptical, whereas the phase-based estimation method is limited to objects that have at least two axes of symmetry. A method based on the extinction theorem was described to estimate the average absorption in the object and demonstrated to give an accurate estimate of absorption coefficient with measured scattering data. For the measured scattering data considered in this investigation, the curve-fitting method gave a more accurate and stable result for the offset in the data than the phase-based method, but phase-based correction of offset performed well in all cases. Each of described methods can be used in appropriate situations to obtain good initial estimates for iterative methods of inverse scattering.

ACKNOWLEDGMENTS

J. C. Tillett is thanked for his suggestions, helpful comments, and useful discussions throughout the duration of this study as well as for performing the scattering calculations that used the k-space method. C. G. Draeger and D. P. Duncan are thanked for the development of codes that used analytic expressions to calculate scattering. A. J. Hesford is thanked for constructive suggestions that improved parts of the narrative. M. M. Foster is thanked for help in the preparation of figures.

The participation of G. Salahura and J. D. Wilbur in reviews of this work as it progressed is acknowledged with appreciation.

This research was funded by NIH Grants EB 000280 and EB 009692 and the University of Rochester Diagnostic Ultrasound research Laboratory Industrial Associates.

Biographies

An external file that holds a picture, illustration, etc.
Object name is nihms214914b1.gif

Jing Jin obtained the B.S. and M.S. degrees from the Department of Applied Physics, South China University of Technology, Guangzhou, China, in 2003 and 2005, respectively. Since September 2005, she has been with the Ultrasound Research Laboratory of the University of Rochester, where she currently is a Ph.D. student in the Department of Electrical and Computer Engineering. Her research interests include image registration, inverse scattering, estimation of scattering object characteristics, and image reconstruction.

An external file that holds a picture, illustration, etc.
Object name is nihms214914b2.gif

Jeffrey P. Astheimer received a B.A. degree in mathematics in 1975, an M.A. degree in mathematics in 1977, and a Ph.D. Degree in mathematics in 1982, all from the University of Rochester. After completing his Ph.D. degree, he joined the faculty of colgate University in 1984, but left in 1985 to co-found the adaptable laboratory software company which developed the Asyst software system. Dr. Astheimer spent 19 years in commercial scientific software development, but has also retained a long-term interest in the mathematics of wave propagation, especially as it relates to applications in medical ultrasound.

An external file that holds a picture, illustration, etc.
Object name is nihms214914b3.gif

Robert C. Waag received his B.E.E., M.S., and Ph.D. degrees from Cornell University in 1961, 1963, and 1965, respectively. After completing his Ph.D. studies, he became a member of the technical staff at Sandia Laboratories, Albuquerque, NM, and then served as an officer in the United States Air Force from 1966 to 1969 at the Rome Air Development Center, Griffiss Air Force Base, NY. In 1969, he joined the faculty of the University of Rochester where he is now Arthur Gold Yates Professor in the Department of Electrical Engineering, School of Engineering and Applied Science, and also holds an appointment in the Department of Radiology, School of Medicine and Dentistry. Prof. Waag’s recent research has treated ultrasonic scattering, propagation, and imaging in medical applications. In 1992, he received the Joseph H. Holmes Pioneer award from the American Institute of Ultrasound in Medicine. He is a life fellow of the IEEE and a fellow of the Acoustical Society of America and the American Institute of Ultrasound in Medicine.

Footnotes

1Because very small distances such as the radius listed in Table I may be accurately estimated in computer simulations to a tenth of a micron, all distances have been reported with this precision.

2For a history of the optical theorem, see [13].

REFERENCES

1. Brignone M, Piana M. The use of constraints for solving inverse scattering problems: Physical optics and the linear sampling method. Inverse Probl. 2005;vol. 21(no. 1):207–222.
2. Kirsch A. Characterization of the shape of a scattering obstacle using the spectral data of the far-field operator. Inverse Probl. 1998;vol. 14(no. 6):1489–1512.
3. Waag RC, Fedewa RJ. A ring transducer system for medical ultrasound research. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2006;vol. 53(no. 10):1707–1718. [PubMed]
4. Lin F, Nachman AI, Waag RC. Quantitative imaging using a time-domain eigenfunction method. J. Acoust. Soc. Am. 2000;vol. 108(no. 3, pt. 1):899–912. [PubMed]
5. Waag RC, Lin F, Varslot TK, Astheimer JP. An eigen-function method for reconstruction of large-scale and high-contrast objects. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2007;vol. 54(no. 7):1316–1332. [PubMed]
6. Mast TD, Nachman AI, Waag RC. Focusing and imaging using eigenfunctions of the scattering operator. J. Acoust. Soc. Am. 1997;vol. 102(no. 2, pt. 1):715–725. [PubMed]
7. Litman A, Lesselier B, Santosa F. Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set. Inverse Probl. 1998;vol. 14(no. 3):685–706.
8. Feijoo GR, Oberai AA, Pinsky PM. An application of shape optimization in the solution of inverse acoustic scattering problems. Inverse Probl. 2004;vol. 20(no. 1):199–228.
9. Hesford AJ, Astheimer JP, Waag RC. Acoustic scattering by arbitrary distributions of disjoint, homogeneous cylinders or spheres. J. Acoust. Soc. Am. 2010;vol. 127(no. 2):850–861. [PubMed]
10. Tabei M, Mast TD, Waag RC. A k-space method for coupled first-order acoustic propagation equations. J. Acoust. Soc. Am. 2002;vol. 111(no. 1, pt. 1):53–63. [PubMed]
11. Lowan A, Morse P, Feshbach H, Lax M. Office of Research and Inventions, Report 62.1R. U.S. Navy Dept.; 1945. Scattering and radiation from circular cylinders and spheres–tables of amplitudes and phase angles.
12. Morse PM, Ingard KU. Theoretical Acoustics. New York, NY: McGraw-Hill; 1968. ch. 8.
13. Newton RG. The optical theorem and beyond. Am. J. Phys. 1976;vol. 44(no. 7):639–642.
14. van de Hulst HC. Light Scattering by Small Particles. New York, NY: Wiley; 1957. p. 39.
15. Born M, Wolf E. Principles of Optics: Electromagnetic Theory of Propagation, Interference, and Diffraction of Light. 7th ed. Cambridge, UK: Cambridge University Press; 2000. p. 720.