акустика / gan_ws_acoustical_imaging_techniques_and_applications_for_en

4

Common Methodologies of

Acoustical Imaging

4.1Introduction

Methodologies of acoustical imaging can be classified into the following categories: (a) tomography, (b) holography, (c) pulse–echo and transmission modes and (d) acoustic microscopy.

4.2Tomography

The concept of acoustical tomography originated from X-ray computer-aided tomography, which is the first successful commercial application of tomography. Tomography is a form of imaging modalities which enables the viewing of the cross-sectional images of an object at various depths, and the combination of these two-dimensional images into a three-dimensional image. The analytical foundation for this method was established in 1917 [1] by Johann Radon, an Austrian mathematician concerned with gravitational theory. The first practical tomographic reconstructions were done 40 years later by Ronald Bracewell [2], a radio astronomer. Radon’s theorem proved that any two-dimensional object can be reconstructed uniquely from an infinite set of its projections. This result has been independently rediscovered a number of times since then by other mathematicians, radio astronomers, electron microscopists, workers in optics, and medical radiologists. The first practical, clinically orientated solutions to reconstructive tomography appeared in the early 1960s after Alan Cormack [3] of Tufts University began to popularize and extend Radon’s work. David Kuhl [4] and co-workers at the University of Pennsylvania built a transverse section scanner for application in nuclear medicine. Kuhl’s scanner was the first tomographic device to isolate for image reconstruction, a single plane transverse to the long axis of the patient’s body. It completely eliminated information from other planes. But the real breakthrough came in 1971 when EMI Ltd in England announced the development of the EMI scanner – a system invented by Godfrey Hounsfield [5] that combined X-ray scanning with digital computing. The system generates images of isolated slices of the

Acoustical Imaging: Techniques and Applications for Engineers, First Edition. Woon Siong Gan. © 2012 John Wiley & Sons, Ltd. Published 2012 by John Wiley & Sons, Ltd.

brain with excellent contrast even for tissues with very small differences in their ability to absorb X-rays.

Tomography was introduced into ultrasonics as an eventual development. The first reported experimental tomograms with ultrasound was due to Greenleaf et al. [6]. It approximates sound wave propagation as straight paths in solids without diffraction. It is much more difficult to take account of diffraction and obtain quantitatively accurate tomograms using ultrasound than to do the same with an X-ray which travels in straight lines in solids with no diffraction. Greenleaf et al. demonstrated the importance of choosing the correct parameter for reconstruction. They showed that, whereas ultrasonic tomography based on absorption shadows was possible, there are inherently many inaccuracies due to reflection, refraction and diffraction. In their method, cross-sections of the object are approximated by N × N grids, with the absorption assumed to be an unknown constant in each block of the grid. A number of projections are then used to solve simultaneous linear equations in the N2 absorption values in each block.

As a remedy to the reflection and refraction problems in absorption-data reconstruction, the reconstruction of the distribution of the ultrasound velocity in the inhomogeneity was proposed [7].

The first papers to deal with diffraction tomography were by Mueller et al. [8]. Nevertheless the properties of ultrasound give rise to several techniques not possible with X-ray. For instance, ultrasound can be subjected to reflection and therefore reflection modes besides transmission modes are also possible with radiation. X-rays on the other hand are hard to reflect and transmission modes of imaging are the only practical modes with their use. Ultrasound is subjected to refraction and travels relatively slowly. Hence the travel time for ultrasound can be easily measured and tomograms of velocity variations as well as attenuation variations can be recorded. Besides this, ultrasound has a coherence property that enables complex amplitudes to be recorded. In order to provide a more comprehensive view and approach to acoustical tomography in this book, we will deal only with diffraction tomography. We will follow the approach of Mueller et al. [8].

For a more accurate representation of the physical processes involved in ultrasonic tomography we shall follow the treatment involving diffraction that takes place when the sound wavelength is much longer than the size of the scatterer. It is only when the sound wavelength is much less than the size of the scatterers that the geometrical optics arguments are valid.

The geometrical optics treatment is also known as a ray-tracing method, which is suitable for X-ray propagation or straight-line propagation where there is no diffraction. When there is a diffraction, the wavefield representation has to be used. In this technique, the acoustic equation of motion or acoustic field equation is reduced to the acoustic wave equation. Here we are dealing only with infinitesimal amplitude sound waves and only the linear acoustic wave equation will be used. This wave equation governs the propagation of sound through a solid medium with spatially varying material parameters. We also take the approximation that the Galilean transformation of the sound propagating system is ignored and, therefore, the time vector is not involved. Even so, only an approximate solution for this wave equation is possible, not an exact analytical solution. The solution is then used to develop an algorithm to reconstruct one or more of the parameters from the observed two-dimensional complex amplitude distribution of sound waves propagated through the medium of interest.

The geometry of the system for ultrasonic tomography using a straight path approximation and the geometric measurement with plane wave insonification are shown in Figures 4.1 and 4.2, respectively.

and F (r ) = 0 for |r| > 0, where F (r ) is the forcing function. The forcing function (4.6) is only valid provided the firstand higher-order derivatives of the medium parameters can be ignored. We now rewrite F (r ) = O(r )ϕ(r) where O(r ) is the object function. O(r ) will be used to represent all inhomogeneities of the object. The object will be reconstructed in terms of the object function O(r ).

We consider the field ϕ(r ) to be the sum of two components. The incident field is the field present without any inhomogeneities or a solution to the homogeneous Helmholtz wave equation:

This leaves the scattered field ϕs (r ) as that part of the field due to the object inhomogeneities, or

The wave equation becomes the inhomogeneous Helmholtz wave equation:

Equation (4.10) is a nonlinear wave equation and cannot be solved for ϕs (r) directly.

In order to solve this equation, it must be linearized. This is normally done in one of two ways: (1) using the Born approximation, (2) using the Rytov approximation.

Prior to this, it is first necessary to write the solution in terms of the Green’s function [11], which is a solution of the differential equation

In two dimensions the solution of (4.11) is written in terms of a zero-order Hankel function of the first kind, and can be expressed as:

The object function in (4.11) represents a point inhomogeneity and the Green’s function can be considered to represent the field resulting from a single point scatterer. Hence it is possible to represent the forcing function of (4.11) as an array of impulses or

O(r )ϕ(r ) = O(r )ϕ(r )δ(r − r )dr (4.14)

The Green’s function represents the solution of the wave equation for a single delta function because the left-hand side of the wave equation is linear and a solution can be written by summing the scattered field due to each individual point scatterer.

Using this idea, the total field due to the impulse O(r )ϕ(r )δ(r − r ) is written as a summation of scaled and shifted versions of the impulse response g(r). This is a simple

convolution and the total radiation from all sources on the right-hand side of (4.10) must be given by the following superposition:

ϕs (r ) = g(r − r )O(r )ϕ(r )dr (4.15)

This is a nonlinear integral equation representing the scattered field ϕs (r ) in terms of the object function O(r ). It is known as a Fredholm equation of the second kind and is a form of scattering integral.

4.2.1The Born Approximation

The Born approximation is derived from high energy physics’ multiple scattering. It is the simpler of the two approximations and is a perturbation approach. The total field is expressed as the sum of the incident field ϕ0 (r ) and a small perturbation, ϕs (r ):

If the scattered field us (r ) is small compared to ϕ0 (r ), the second integral in (4.17) can be ignored and

4.2.2The Rytov Approximation

The Rytov approximation is another first-order approximation to the scattered field and is valid under different restrictions. It is also a perturbation method.

Consider the total field to be represented as a complex phase [12],

Expressing the total complex phase ψ as the sum of the incident phase function ψ0 and the scattered complex phase, ψs (r ), we have

Consider the first Rytov approximation to the scattered field, then

Substituting in ϕ0 = Ae jk0s·r for the incident field and in ψs = exp( − jk0s · r)ϕB (r ) for the scattered phase, then

For small ϕB, expanding the first exponential in terms of its power series and keeping only the first two terms, we have

Thus, for very small objects and perturbation, the Rytov solution is approximately equal to the Born solution given in (4.18).

The first-order Born and Rytov solutions will provide information on the scattered field to be used in the reconstruction. In the Born approximation, the complex amplitude of the scattered field is measured and used as an estimate of the function uB, while in the Rytov approximation uB is estimated from the phase of the scattered field. The Rytov approximation is considered to be more accurate than the Born approximation and should provide a better estimate of uB.

4.2.3The Fourier Diffraction Theorem

This theorem [13], which is fundamental to diffraction tomography, relates the Fourier transform of the measured forward-scattered data with the Fourier transform of the object. The theorem is valid when the inhomogeneities in the object are only weakly scattering. The theorem states that when an object O(x, y) is illuminated with a plane wave, as shown in Figure 4.3, the Fourier transform of the forward-scattered fields measured on the line TT gives the values of the two-dimensional Fourier, O(ω1,ω2) of the object along a circular arc in the frequency domain, as shown in the right half of the figure.

The importance of this theorem is made obvious by noting that if an object is illuminated by plane waves from many directions over 360◦, the resulting circular arcs in the (ω1,ω2) plane will fill up the frequency domain. The function O(x, y) can then be reconstructed by an inverse Fourier transform. This theorem relates to any cylindrical object whose cross-sectional distribution is given by the function O(x, y). The forward-scattered fields are measured on a line

Figure 4.3 The Fourier diffraction theorem (Kak and Slaney [14] © IEEE)

of detection along TT in Figure 4.3. If a truly three-dimensional object were illuminated by the plane wave, the forward-scattered fields would now have to be measured by a plane array of detectors. The Fourier transform of the fields measured by such an array would give the values of the three-dimensional transform of the object over a spherical surface [15]. Subsequently, works by Nahamoo et al. [16] and Devaney [17] showed a new synthetic aperture procedure for a full three-dimensional reconstruction using only two rotational positions of the object.

4.2.4Reconstruction and Backpropagation Algorithm

The Fourier Diffraction Theorem shows that when an object is illuminated with a plane wave travelling in the positive y direction, the Fourier transform of the forward-scattered fields gives

values of the object’s Fourier transform on an arc. Therefore, if an object is illuminated from

√

many different directions it is possible, in principle, to fill up a disc of diameter 2k in the frequency domain with samples of the Fourier transform of the object, and then reconstruct the object by direct Fourier inversion [18]. Therefore, diffraction tomography, using forward-

scattered data only, determines the object up to a maximum angular spatial frequency of

√

2k. To this extent, the reconstructed object is a low-pass version of the original. In practice, the loss of resolution caused by this bandlimiting is negligible, being more influenced by considerations such as the aperture sizes of the transmitting and receiving elements.

The fact that the frequency domain samples are available over circular arcs (whereas for convenient display it is desired to have samples over a rectangular lattice) is a source of

computational difficulty in reconstruction algorithms for diffraction tomography. It should also be pointed out that by illuminating the object over 360◦, a uniform double coverage of the frequency domain is generated. If the illumination is restricted to a portion of 360◦ there will still be a complete coverage of the frequency domain, but in this case there will be patches in the (ω1,ω2) plane that would have a double coverage. In reconstructing from circular arc grids to rectangular grids, it is often easier to contend with a uniform double coverage rather than a coverage that is single in most cases and only double in patches.

However, for some applications not given to data collection from all possible directions, it is useful to bear in mind that it is not necessary to go completely around an object to get complete coverage of the frequency domain. In principle, it should be possible to get an equal quality reconstruction when illumination angles are restricted to the 180◦ required to complete the coverage of the frequency domain.

There are two computational strategies for reconstructing the object given measurements of the scattered field. Soumekh et al. [19] pointed out that the two algorithms can be considered as interpolation in the frequency domain and in the space domain, and are analogous to the direct Fourier inversion and backprojection of X-ray tomography. But unlike X-ray tomography, where backprojection is the preferred approach, the computational expense of a space domain interpolation of diffracted projections makes the frequency domain interpolation the preferred approach. We will focus only on the frequency domain interpolation. First we will discuss the frequency domain interpolation between a circular grid in which the data is generated by diffraction tomography, and a rectangular grid suitable for image reconstruction, where the parameters representing each grid must be selected and the relationship between the two sets of parameters must then be written.

The backpropgation algorithm for reconstruction was introduced by Devaney [20] and Kaveh [21]. This is analogous to the filtered-backprojection algorithm used in X-ray tomography. The filtered-backprojection algorithm contributes to the success of X-ray tomography due to its superior numerical accuracy. Unfortunately, a backpropagation algorithm, on the other hand, does not possess the same efficiency in implementation as it is much more computationally intensive. With regard to accuracy, it does not seem to possess any advantage, especially if interpolation is carried out after increasing the sampling density by appropriate zero-padding.

The derivation of the backpropagation algorithm will follow the procedure of Devaney [20]. The inverse Fourier transform of the object function is given by

αα

This integral represents the object function in terms of the Fourier transform of the object along a rectangular grid. Since a diffraction tomography experiment measures the Fourier transform of the object along circular arcs, it will be easier to perform the integration if it is modified to use the projection data more naturally. This can be done using two coordinate transformations: the first will exchange the rectangular grid for a set of semicircular arcs, and the second will map the arcs into their plane wave decomposition.

To exchange the rectangular grid for semicircular arcs, we represent