## Abstract

Tomography of complex three-dimensional objects with ultrasound or microwave has been a long-standing goal since the introduction of these technologies after World War II. While current state-of-the-art systems can provide high-resolution images of cylindrical objects characterized by a two-dimensional structure, the three-dimensional case remains an open challenge owing to current limitations of sensor technology and computer power. Here, this problem is addressed by means of a synthetic aperture technique that, while using hardware technology developed for two-dimensional problems, accounts for the complexity of three-dimensional scattering and leads to high-resolution three-dimensional reconstructions. In this paper, we present the theoretical formulation of this new approach and illustrate it by means of a numerical example.

## 1. Introduction

Tomography attempts to reconstruct the spatial distribution of one or more physical parameters of an object by studying the perturbation induced by the object’s structure on the free propagation of either mechanical or electromagnetic waves. The wave–matter interaction can be modelled according to classical ray theory or under the more general framework of diffraction theory, which includes the former in the short wavelength limit. Thus, while a ray is, in general, sufficient to describe the propagation of high-energy photons in X-ray tomography of biological materials, diffraction, refraction and multiple scattering can become dominant when imaging the same material with ultrasound or microwave. The presence of these effects poses a number of fundamental challenges to the development of diffraction tomography (DT) technology. To appreciate the complexity of the problem, we now consider the application of ultrasound tomography to breast cancer detection that has been investigated since the early 1970s (Greenleaf *et al.*1973, 1975; Carson *et al.* 1981) and is yet to be used in clinical diagnostics.

Breast is arguably one of the most complex tissues in the human body, and despite the fact that its anatomy is not completely understood (Ramsay *et al.* 2005), it is largely accepted that the gross anatomy consists of a fatty tissue and a glandular component that is arranged in lobes that are held together by a fibrous framework called the Cooper’s ligaments. The lobes contain lobules that are connected to the nipple by an entangled network of ducts that Cooper, in the celebrated *Anatomy of the breast*, compared with the intertwined roots of a tree; in addition, lymphatic and blood vessels provide further complexity to the architecture of the breast. When an ultrasonic wave travels through these structures, their different mechanical properties (mainly density and stiffness) trigger diffraction and scattering phenomena that change the structure of the wave field. The extent of these changes largely depends on the size of the structures relative to wavelength λ of the probing wave. In conventional medical ultrasound, wavelengths range between 0.1 and 1 mm, which correspond to the scale of a wide spectrum of the breast structural features. As a result, the use of ray theory to describe the interaction of ultrasound with breast tissue is not always valid.

The possibility of developing an ultrasound tomography technology that could lead to images of the breast comparable to those obtained with high-resolution X-ray computed tomography (CT) is very attractive. Indeed, ultrasound is devoid of risk for patients (excessive exposure to ionizing radiation precludes the use of CT for *in vivo* breast imaging) and has higher sensitivity than mammography (currently the standard method of screening) in dense breast (Kolb *et al.* 2002; Berg *et al.* 2008). Moreover, recent discoveries in the field of cell and tissue biomechanics are laying the foundations to improve our current understanding of the links between ultrasound measurements and changes in tissue mechanical properties that mediate the progression of cancer (Wang *et al.* 1993; Huang & Ingber 2005). In this context, the system developed by Andre *et al.* (1997) and, more recently, the computed ultrasound risk evaluation (CURE) prototype (Duric *et al.* 2007) are examples of ultrasound tomography setups for breast imaging. Instead of scanning the breast with a freehand probe as in conventional sonography (Wells 1999), a patient lies prone on a table with a breast suspended in a water bath through an aperture in the table (figure 1). A toroidal array, probe then encircles the breast and is scanned vertically from the chest wall to the nipple region to sweep the entire volume of the breast. At each array position, the breast is insonified from all possible angles in the plane of the array, and the backscattered and transmitted ultrasonic fields through the breast are measured by all the array elements in parallel and for each insonification. These measurements are subsequently used by a tomography algorithm to reconstruct a slice across the breast in the plane of the array. This results in a fully automated scanning procedure that is operator independent and does not require the compression of the breast, which is known to cause apprehension in women undergoing mammography (Kornguth *et al.* 1996).

Different algorithms have been proposed to reconstruct slices across the breast, such as ray-based travel time tomography (Greenleaf *et al.* 1973; Carson *et al.* 1981; Duric *et al.* 2007) and DT (Devaney 1983). These techniques are linear inverse methods that provide a filtered reconstruction of the true image. Therefore, inverse filtering methods (Bertero & Bocacci 1998) have been considered for image restoration such as the filtered back-propagation method for DT introduced by Devaney (1982). Current applications of tomography to breast imaging have been based on two-dimensional reconstruction methods that image a slice from the measurements performed by the array at a particular position and assume that the ultrasonic field ‘sees’ the slice corresponding to the middle plane of the array only. The use of the two-dimensional approximation rather than full three-dimensional models has been dictated by ultrasonic hardware limitations and high computational cost (Solimene *et al.* in press). While the two-dimensional approximation is realistic in the case of X-ray CT as photons travel along straight trajectories, in the case of ultrasonic probing it is no longer possible to confine the region of interaction between the wave and the object to a straight ray because of diffraction. This implies that the field measured along the array aperture is not only dependent on the properties of the object within the plane of the aperture, but is also affected by structures outside it. As a result, image artefact and resolution degradation occur (Williamson & Worthington 1993), particularly in the direction perpendicular to the plane of the array. The poor resolution in the scanning direction increases what is known as the slice thickness and severely limits the three-dimensional imaging capability of ultrasound tomography.

Undoubtedly, breast tomography remains the major potential application of ultrasound tomography; however, the issues discussed earlier are common to other applications of DT in optics, geoscience and non-destructive testing when objects with complex internal structures are probed with mechanical or electromagnetic waves. In these cases, the wavelength required for the ray approximation to be valid must be very short compared with the internal structures of the object. However, the same structures cause significant scattering outside the plane of the array that attenuates the strength of the signal propagating across the object dramatically, leading to weak signals that are difficult to detect. Moreover, the same structures can destroy the coherence of the wave field, which would then propagate as a diffusive field (Page *et al.* 1995; Weaver & Lobkis 2006).

This paper introduces a new imaging method for three-dimensional DT that uses the same measurement setup depicted in figure 1 while considering three-dimensional scattering instead of the two-dimensional models used in previous papers. This method is inspired by synthetic aperture radar techniques used in spaceborne imaging of the Earth’s surface. Here, high-resolution images of the Earth’s topography are obtained by combining the signals transmitted by a microwave antenna and reflected by the Earth’s surface as the antenna moves along the trajectory of a satellite (Bamler & Hartl 1998). This leads to a synthetic aperture that is greater than the actual aperture of the antenna and improves the resolution along the direction of the antenna path. In this paper, the same concept is used to improve the resolution in the array scanning direction and leads to the concept of three-dimensional synthetic aperture DT (SADT). This term was first introduced by Nahamoo *et al.* (1984) in a different context of two-dimensional DT, where two parallel linear array transducers were mimicked by means of two single element transducers scanned along straight lines.

The paper begins by reviewing the classical scattering models used in DT to link the measurements and the structure of the object represented in the spatial frequency domain in §2. Next, §3 introduces the ‘ideal’ approach to three-dimensional DT and discusses the reasons why this is not feasible in practice. Section 4 presents the new SADT method that is subsequently validated by means of a numerical example in §5. Section 6 discusses practical experimental issues for the implementation of SADT and is followed by conclusions in §7.

## 2. Measurements and object’s structure

Let us consider a complex three-dimensional object that occupies a volume *D* in the space denoted by a Cartesian coordinate system {*O*,*x*_{1},*x*_{2},*x*_{3}} as shown in figure 2*a*. In this paper, it is assumed that the scattering problem is described by a scalar wave field, *ψ*, solution to
2.1
where is the Helmholtz operator (), *k*_{0} is the background wavenumber (2* λ*/λ), specifies the direction of an incident plane wave that illuminates the object and

*ω*is the angular frequency. The unit vector is defined by the angles

*θ*

_{t}and

*ϕ*

_{t}. The object is described by the so-called

*Object Function*that depends on the type of wave field used to probe the object: for electromagnetic wave sensing, it is related to the index of refraction (Born & Wolf 1999),

*n*(

**r**,

*ω*), through the relation , and for acoustic waves, it is linked to the speed of sound and the attenuation coefficient (). In particular, for a lossless object 2.2 where

*c*

_{0}is the sound speed of the homogeneous background in which the object is immersed and

*c*(

**r**,

*ω*) is the local sound speed inside the object. The dependence of the object function on

*ω*is because of dispersion and energy dissipation phenomena. The analysis performed in the rest of this paper will consider monochromatic wave fields; therefore, the explicit dependence on

*ω*is omitted.

The DT problem consists of reconstructing the function *O*(**r**) from a set of scattering experiments. For this purpose, it is convenient to introduce the representation of the object function in the spatial frequency domain, K-space, which is obtained by performing the three-dimensional Fourier transform of *O*(**r**)
2.3

It is now observed that the total field far away from the object can be expressed as (Born & Wolf 1999) 2.4 where is the scattering amplitude defined as 2.5 Through equation (2.4), the scattering amplitude can be measured experimentally, by illuminating the object from all possible directions and for each of them detecting the scattered field in all directions , both illumination and detection being performed in the far field.

Under the Born approximation, the fundamental Fourier diffraction theorem (Wolf 1969) demonstrates that the scattering amplitude corresponding to a particular illumination and detection pair, , provides at the spatial frequency , i.e.
2.6
as shown in the diagram of figure 2*b*. Therefore, if scattering measurements could be performed for all possible pairs and , could be reconstructed for all the spatial frequencies within a ball of radius 2*k*_{0}, also known as the Ewald limiting sphere (ELS) (Wolf 1969). A similar mapping is obtained under the Rytov approximation, in which instead of the scattering amplitude in equation (2.6), the so-called Rytov data are used (Sponheim *et al.* 1991).

The linear mapping is central to DT algorithms and is a direct consequence of the linear approximations introduced by the Born and Rytov hypotheses. Indeed, as shown in Simonetti (2006), the mapping becomes nonlinear when more complex scattering models that include multiple scattering are considered. Moreover, it can be observed that the linear approximation implies that some scattering measurements are redundant. For instance, all the through-transmission measurements () map onto the origin of the K-space, meaning that through-transmission measurements are the same regardless of the illumination direction. Obviously, this is not physical when the object lacks spherical symmetry and leads to a paradox. We will come back to this issue in §4 when SADT is introduced.

The linear mapping in equation (2.6) provides a means to reconstruct *O*(**r**). Indeed, *O*(**r**) can be retrieved by performing the inverse Fourier transform of the reconstructed , assuming that it vanishes outside the ELS. This leads to a low-pass-filtered reconstruction of the object function containing spatial frequencies smaller than 2*k*_{0} according to the classical resolution limit, λ/2. Although this approach is conceptually straightforward, its practical implementation for three-dimensional imaging is very challenging. The main difficulty lies in the number of scattering experiments required to populate the ELS in order to achieve sufficient resolution and avoid aliasing artefact.

## 3. The ideal spherical aperture

In this section, the feasibility of three-dimensional DT is discussed after introducing an alternative approach to standard DT algorithms such as the filtered backpropagation method (Devaney 1982), based on a three-dimensional beamforming (BF) algorithm. This algorithm also provides the basis for the SADT technique developed in §4.

### (a) Three-dimensional beamforming algorithm

Let us assume that the scattering amplitude, *f*(*θ*_{r},*θ*_{t},*ϕ*_{r},*ϕ*_{t}), can be measured as a continuous function of the illumination and detection directions, i.e. *θ*_{r},*θ*_{t}∈[0,2* λ*] and

*ϕ*

_{r},

*ϕ*

_{t}∈[0,

*], these angles being shown in figure 2*

*λ**a*. In principle, this could be achieved with a spherical array of transreceivers that surrounds the object.

Standard BF produces the image of an object at a point, **z**, of the image space by focusing an incident beam at **r**=**z** in the object space. The resulting scattered field is subsequently phase shifted and integrated over the aperture of the array, so that only the contributions to the scattered field from the focal point are added coherently. This two-step process is obtained by means of the BF functional
3.1
where is the unit vector associated with the angles *θ* and *ϕ*. As discussed by Simonetti & Huang (2008) for the two-dimensional case, the second exponential in equation (3.1) represents focusing in transmission, whereas the first corresponds to the focusing of the received scattered field. The point spread function (PSF) associated with the functional (3.1) can be obtained by considering the image of a point scatterer at position **r**. In this case, the scattering amplitude is
3.2
and the PSF reads
3.3
which, with a straightforward calculation, leads to
3.4
where *j*_{0} is the spherical Bessel function of the order 0. The dependence of *h*_{BF} on |**z**−**r**| rather than on the individual variables, **z** and **r**, implies that the PSF is space invariant. The Fourier transform of *h*_{BF} becomes
3.5
Under a linear scattering model and using the spatial invariance of *h*_{BF}, the BF image of the object function *O*(**r**) is given by the convolution integral
3.6
that in the spatial frequency domain reads
3.7
While DT leads to the low-pass-filtered image, , the BF algorithm introduces a distortion that is described by the additional filter 1/|**k**|. As a result, the DT image can be obtained from the BF image by applying the filter |**k**| to the BF image. This is an alternative approach to other DT algorithms ().

### (b) Practical challenges

Regardless of the algorithm used to form the image, three-dimensional DT requires a vast amount of data. The key parameters determining the volume of data are the wavelength of the probing field and the characteristic size of the object, here defined as the radius, *R*_{0}, of the sphere that circumscribes the object. These parameters define the sampling condition required to capture the information contained in the scattered fields. This can be shown by considering the scattered field that radiates from the sphere that circumscribes the object to an ideal spherical array of detectors of radius *R* concentric with the sphere, i.e.
3.8
where *a*_{mn} are coefficients that depend on the distribution of the scattered field over the sphere of radius *R*_{0}, is the spherical Hankel function of the first kind and order *n* representing outgoing waves and
3.9
are the spherical harmonics of order *n* and degree *m* that depend on the Legendre polynomials, . Assuming that *R*>>*R*_{0}, only the orders give a significant contribution to the expansion in equation (3.8) owing to the asymptotic behaviour of the spherical Hankel functions. As a result, the scattered field can be considered as bandlimited with bandwidth . The sampling criterion for such a field was given by Driscoll & Healy (1994) and states that the full information can be retrieved with *B*^{2} sampling points distributed over an equiangular grid of points (*ϕ*_{i},*θ*_{j}), *i*,*j*=0,…,2*B*−1, where *ϕ*_{i}=*λ**i*/2*B* and *θ*_{j}=*λ**j*/*B*. Based on a reciprocity argument, the illumination directions also have to cover the same equiangular grid. In principle, this could be achieved by deploying ultrasonic transducers, which can act as both transmitters and receivers, along the spherical aperture. However, in practical applications, the number of transducers is usually unrealistic. As an example, for an object with a characteristic size-to-wavelength ratio of the order of 100 (this is realistic in the case of breast imaging), the number of sampling points would be in the region of approximately equal to 400 000. This is beyond the capabilities of current ultrasound technology because of limitations in array manufacturing and front-end electronics. Moreover, the number of scattering experiments to be considered is of the order of *B*^{4}≈10^{10}, which is unfeasible owing to the amount of memory and computer power needed to process the data as well as the time required to perform such a vast amount of measurements. To overcome these limitations, Halse *et al.* (2003) have proposed to illuminate the object from a set of directions perpendicular to a prescribed axis and to measure for each of them the transmitted field over a planar array of receivers. However, the reconstruction suffered from distortion artefacts, particularly in those portions of the object distant from the plane perpendicular to the axis and passing through the centre of the array. Moreover, such a method discards all the information contained in backscattered waves.

## 4. Synthetic aperture diffraction tomography

For the reasons discussed in §3, current applications of ultrasound tomography make use of the ray approximation that allows the object to be imaged slice by slice with a toroidal array (figure 1). As mentioned previously, the ray approximation causes artefact and leads to image resolution that is poorer than the theoretical λ/2 limit. More sophisticated models that account for diffraction have only been investigated for cylindrical objects that can be treated as two dimensional (Lin *et al.* 2000; Simonetti *et al.* 2007; Waag *et al.* 2007). The absence of the third dimension eliminates one degree of freedom in the measurements, which can then be performed with a toroidal array rather than a spherical one. A toroidal array architecture can now be built owing to recent progress in solid-state electronics and array micro-machining that allows a large number of scattering experiments to be performed in a very short period of time. For instance, prototype toroidal arrays with as many as 1024 (Andre *et al.* 1997), or more recently, 2048 (Waag & Fedewa 2006), transreceivers have been manufactured and can collect all the possible transmit–receive pairs in less than a second. Acquisition speed is crucial for clinical applications to avoid motion artefact and limit examination time.

The SADT method proposed in this paper employs this same technology; however, a slice across the object is now reconstructed by combining the measurements obtained for all the array positions rather than the single dataset measured when the array intersects the slice. Before introducing the method, we show that under a linear scattering model, the measurements performed by scanning a toroidal array are sufficient to populate the ELS.

With reference to figure 3*a*, let the elevation angle *ϕ* denote the vertical position of the toroidal array. The illumination and detection directions are now identified by the angle pairs (*θ*_{t}, *ϕ*) and (*θ*_{r}, *ϕ*), respectively. At a particular array position *ϕ*, the components of the spatial frequency vector, , are
4.1
4.2
and
4.3
therefore, as the (*θ*_{t},*θ*_{r}) pairs vary in the interval [0,2* λ*]×[0,2

*], the measurements populate a disc of the K-space, contained in the plane and with radius as shown in figure 3*

*λ**b*. As a result, by scanning the array vertically, the entire ELS can be filled despite the fact that not all the and pairs are available. This is consistent with the observation that under a linear scattering model, many of the possible scattering experiments are redundant. Therefore, a full-view configuration is not necessary to retrieve all the information, and the measurements performed in a synthetic aperture fashion are sufficient to reconstruct the object function as in the ideal case of a spherical array of sensors. Importantly, spatial frequencies as high as 2

*k*

_{0}can be retrieved in the scanning direction, leading to an axial resolution limit of λ/2.

The SADT algorithm proposed in this paper consists of two steps. First, a synthetic aperture BF (SABF) image of the object is produced. Subsequently, the DT image is obtained from the SABF image by deconvolution with the PSF of SABF. For this purpose, we introduced the SABF algorithm first.

### (a) Synthetic aperture beamforming

To simplify the mathematical treatment, it is assumed that the scattering amplitude *f*(*θ*_{r},*θ*_{t},*ϕ*) can be measured for any *ϕ*∈[0,* λ*], i.e. the array is scanned from to . In addition, for each position of the array, all the transmit–receive pairs are measured, i.e.

*θ*

_{r},

*θ*

_{t}∈[0,2

*]. Let*

*λ***z**be a point in the image space, we define the BF functional in analogy with equation (3.1) as 4.4 where only one integral over the elevation angle is performed. The PSF associated with the functional (4.4) can be evaluated by considering the scattering amplitude of a point scatterer (3.2) and solving the integral 4.5 This is space invariant since it depends on

**z**−

**r**rather than on the individual variables

**z**and

**r**. As a result, the SABF image is the convolution of the object function with the PSF. To evaluate the PSF, let

**a**=

**z**−

**r**. Expression (4.5) can be rewritten in the form 4.6 where the subscripts ∥ and ⊥ refer to the components of a vector in the plane

*x*

_{3}=0 and along , respectively. By using the Jacobi–Anger expansion, the integral in equation (4.6) becomes 4.7 where

*J*

_{0}is the Bessel function of the first kind of zero order. Since

*J*

_{0}is even and is odd, the integral in equation (4.7) is equivalent to 4.8 with α=2

*k*

_{0}

*a*

_{⊥}and β=

*k*

*a*

_{∥}. Equation (4.8) implies that the PSF is real. By using the series expansion of and the Bessel second integral (Watson 1952), the integral in equation (4.8) can be evaluated as a series of Bessel functions of the first kind,

*J*

_{m}, i.e. 4.9 where (2

*m*

*m*) are the binomial coefficients.

The three-dimensional structure of the normalized PSF is shown in figure 4. As can be deduced from equation (4.9), the PSF is axisymmetric relative to the -axis (figure 4*d*) and is described by the function 4*λ*^{2}*J*_{0}(2*k**x*_{3}) along it (figure 4*c*). If the resolution in the direction is defined by considering the first null of the PSF, according to the Rayleigh criterion, the resolution limit becomes 0.19λ at the expense of large sidelobes. On the other hand, the sidelobes are lower in the plane *x*_{3}=0; however, the width of the central lobe is now larger (figure 4*e*).

To estimate the resolution of SABF in a more rigorous fashion, its spatial Fourier transform is considered. This can be calculated from the integral in equation (4.8) as described in appendix A and reads
4.10
As a result, the SABF image is
4.11
The presence of the band-pass filter *Π* implies that the SABF reconstruction contains spatial frequencies up to 2*k*_{0} only, and therefore its resolution is constrained by the classical resolution limit λ/2. Moreover, SABF introduces a significant image distortion to the ideal DT image , described by the factor . SABF amplifies the spatial frequencies corresponding to the shell of the ELS where becomes singular, as well as the spatial frequencies on the *k*_{3}-axis, i.e. **k**_{∥}=0.

Expression (4.11) is the main result of this paper and provides the means to reconstruct a DT image from the SABF image according to 4.12

## 5. Numerical implementation of synthetic aperture diffraction tomography

The numerical implementation of the SADT method is illustrated in the flowchart shown in figure 5. After selecting the volume to be imaged, a three-dimensional cubic grid of points is defined. As it is expected that the highest spatial frequency is 2*k*_{0}, the size of each voxel has to be smaller than λ/4, to satisfy the Nyquist sampling criterion. Subsequently, the SABF functional is estimated at each point of the grid using equation (4.4) and the measured dataset, i.e. the sampled *f*(*θ*_{r},*θ*_{t},*ϕ*). The DT image is subsequently obtained by deconvolving *h*_{SABF} from the SABF image. This is most efficiently achieved in the spatial frequency domain by performing the three-dimensional fast Fourier transform (FFT) of the SABF image, dividing it by *H*_{SABF} and performing the inverse FFT (IFFT) to go back to the geometric space. In this process, the most computationally extensive step is the formation of the three-dimensional SABF image. However, it should be observed that as long as the imaging grid is uniform, there is no need for data interpolation techniques as required by standard DT algorithms, thus resulting in faster reconstructions.

A numerical example is shown next. Here, the reconstruction of the speed of sound distribution inside a homogeneous sphere is considered. In this case, the forward scattering can be predicted with a semi-analytic formula (see Anderson 1949), which does not use a linearized forward scattering model and accounts for multiple scattering inside the sphere. The speed of sound inside the sphere is set to be *c*=1485 ms^{−1} and the sphere is immersed in water background (*c*_{0}=1500 ms^{−1}), this level of contrast being typical of human tissue. Moreover, no density contrast is considered. The radius of the sphere is 20 mm and is insonified by plane waves at 100 KHz, i.e. λ=15 mm. The toroidal array has a radius of 200 mm and consists of *N*=80 point-like transducers, and measurements are performed over 79 vertical positions. The vertical positions are selected so as to achieve a constant elevation interval Δ*ϕ*=* λ*/

*N*. The imaging volume is a cube 248 mm in size and divided into 1 mm cubic voxels.

Figure 6 shows the three-dimensional SABF image of the sphere. Due to the anisotropy of the PSF, *h*_{SABF}, the reconstruction of the sphere is not homogeneous as it is clearly seen in figure 6*b*,*d* that are two cross sections in the vertical and horizontal planes passing through the origin of the sphere. The reconstruction along the -axis (figure 6*c*) is sharper than that in the *x*_{3}=0 plane (figure 6*e*); however, it contains more artefacts. The higher resolution is because of the narrow central lobe of the PSF along the -axis (cf. figure 4*c*,*e*). On the other hand, the artefact peaks seen in figure 6*c* are due to the sidelobes of the PSF, which are larger along the -axis than those in the *x*_{3}=0 plane (cf. figure 4*c*,*e*).

Next, figure 7 is the DT reconstruction of the object function obtained after deconvolution of the SABF image with *h*_{SABF}. The reconstruction is almost spherically symmetric, and most of the artefacts of the SABF image along the -axis are suppressed (some residual artefact is due to the fact that the scan does not contain the elevation angles *ϕ*=0 and * λ*/2). Moreover, the image is sharper as can be seen in the cross sections in figure 7

*c*,

*e*and also provides an accurate reconstruction of the sphere shape. Using relationship (2.2), the sound speed map can be obtained from the object function as shown in figure 8, which provides the reconstructed sound speed along the -axis.

## 6. Discussion

Under the linear scattering models considered in this paper, DT provides the ideal reconstruction of the object function as the same weight is given to all the spatial frequencies contained within the ELS, as described by the low-pass filter *Π*(|**k**|). On the other hand, different BF algorithms introduce some degree of image distortion as shown in table 1, which compares the distortion filters of BF and SABF. The distortion filter for the two-dimensional problem (Simonetti & Huang 2008) is also provided for comparison. While the distortion filter is isotropic for three-dimensional BF, the filter associated with SABF leads to a strong anisotropy along the *k*_{3}-axis owing to the **k**_{∥} term. The two-dimensional BF algorithm corresponds to the two-dimensional version of the full-view three-dimensional BF discussed in §3. As for three-dimensional BF, the distortion filter is symmetric. However, the structure of the two-dimensional BF filter is closer to that of SABF; in particular, it is the same for the spatial frequencies in the plane *k*_{3}=0 (**k**_{∥}=|**k**|).

The practical implementation of the SADT method is based on the possibility of measuring the scattering amplitude *f*(*θ*_{r},*θ*_{t},*ϕ*) for as wide a range of elevation angles as possible (ideally *ϕ*∈[0,* λ*]). It has to be observed that the travel distance, 2

*t*, of the array is limited; therefore, the actual range of elevation angles is reduced by an amount dependent on the radius of the array,

*R*, so that

*ϕ*∈[ϵ,

*−ϵ]. This causes the spatial frequencies with to be lost and leads to a partial resolution degradation along the scan direction.*

*λ*Another important aspect for the feasibility of SADT is that the solid angle of the illumination beam from each array element should be as wide as possible, approaching the field produced by a point source. If this condition is met, each array element will also be equally sensitive to wavefields incoming from different directions by reciprocity. These two characteristics are crucial to illuminate planes of the object that are distant from the plane of the array and to detect the field scattered from the structures in those planes. A second key parameter is the radius of the array that should be much larger than the object characteristic size. This condition is required to be consistent with the far-field approximation and to ensure that incident beams can be considered as plane waves. The latter requirement can be achieved by ensuring that *R*>*D*^{2}/8λ, where *D* is the characteristic size of the object. This condition guarantees that the phase delay accumulated along the curved wavefront of the incident field and relative to an ideal plane wavefront is less than λ across the object.

To approximate the point-source radiation pattern, the active area of each array element should be a disc with subwavelength diameter. However, radiation efficiency from subwavelength apertures is low, resulting in poor signal-to-noise ratios. This is further reduced by the very larger electro-mechanical impedance of small transducers that causes severe mismatch with the impedance of the electronics used to control the array (typically 50 Ω). This suggests that the use of small aperture elements would not be feasible in practice. An alternative is to use larger array elements, each with a spherical diverging lens that broadens the beam from each element. Indeed, as the number of transducers is dictated by the size of the object relative to the wavelength only, large array elements could be accommodated in a ring of sufficiently large radius.

Finally, with respect to breast cancer detection, the patient chest wall will be the main issue. With reference to figure 3, it can be assumed that the chest wall lies over the plane *x*_{3}=0 and that the breast occupies the space *x*_{3}<0. This implies that only the elevation angles *ϕ*>* λ*/2 can be measured. Moreover, it is necessary to consider appropriate boundary conditions along the plane

*x*

_{3}=0. One possibility is to use the sound-soft boundary condition that simplifies the problem dramatically. In this case, the scattering problem is equivalent to the scattering scenario obtained by considering the actual breast and array and their mirror image relative to

*x*

_{3}=0 and the continuity of the wave field across

*x*

_{3}=0. This condition also allows the information relative to the elevation angles in the interval [0,

*/2] to be retrieved. While the sound-soft condition is exact along the free surface of water, it is not accurate along that portion of the*

*λ**x*

_{3}=0 plane that intersects the breast, and the validity of this approach will have to be tested experimentally.

## 7. Conclusions

A new method for DT of three-dimensional objects has been introduced. It employs a toroidal array of transreceivers that encircles the object and is scanned so that its middle plane sweeps the entire volume of the object. For each array position, all the possible transmit–receive pairs are measured. It has been shown that these measurements are sufficient to populate the ELS as in the ideal case of a spherical array of transducers that surrounds the object. Therefore, it is possible to reconstruct the object structure with the same degree of accuracy as in full-view DT, leading to the concept of SADT.

An SADT algorithm to reconstruct the object function from the measurements has been proposed. It is based on an SABF method that is also introduced in this paper. Noticeably, SADT achieves the theoretical λ/2 resolution limit in the scanning direction, thus limiting the slice thickness. Moreover, SADT is more attractive than full-view DT as it reduces the number of array elements required to avoid aliasing, hence the complexity of the hardware, and limits the size of the datasets to be processed.

The SADT approach has been validated by means of numerical simulations that have considered a homogeneous fluid sphere immersed in a water background. Accurate shape and sound speed distribution inside the sphere have been reconstructed.

## Acknowledgements

This work was supported through the U.S. DOE Laboratory Directed Research and Development program at Los Alamos National Laboratory. F.S. acknowledges the support of the UK Royal Academy of Engineering/EPSRC.

## Appendix A. Spatial frequency analysis of *h*_{SABF}

The three-dimensional Fourier transform of *h*_{SABF} can be calculated from the integral expression in equation (4.8). For this purpose, it can be observed that the three-dimensional Fourier transform of a function with circular symmetry, *q*(*a*_{∥},*a*⊥)*J*_{0}(**k**_{∥}*a*_{∥}), is symmetric and can be written as
A 1
Moreover, if *q*(*a*_{∥},*a*⊥)=*g*(*a*_{∥})*f*(*a*⊥), then
A 2
where is the Hankel transform
A 3
and is the mono-dimensional Fourier transform
A 4
In order to obtain the Fourier transform of the PSF, we apply expression (A 1) to equation (4.8)
A 5
Moreover,
A 6
and
A 7
Using the sampling property of the delta function δ, the sole angular value that contributes to the integral in equation (A 5) is leading to
A 8

## Footnotes

- Received March 31, 2009.
- Accepted June 2, 2009.

- © 2009 The Royal Society