Acoustic interface contrast imaging

In acoustic reflectivity imaging, we infer the internal reflectivity of an unknown object from reflected waveforms. A common assumption is that the mass density is constant and that the recorded pressure field is related to a volume contrast in the wave speed by a nonlinear volume-integral representation. This representation is typically linearized under the Born approximation and solved for the volume contrast by iterative inversion. We propose an alternative methodology, which we refer to as interface contrast imaging. In our derivation, we assume a medium with constant wave speed, which contains discontinuities of the acoustic impedance at a collection of interfaces between piecewise-homogeneous subdomains. A linear relationship is established between the recorded data and the gradient of the acoustic impedance at the interfaces, which we refer to as an interface contrast. This contrast can be solved for by iterative inversion. With this procedure, acoustic interfaces can be delineated with superior resolution compared to volume contrast imaging. Since the convergence speed is relatively fast and a reasonable image can already be obtained after a single iteration, real-time applications seem feasible. If necessary, the acoustic impedance can also be imaged by integrating the retrieved reflectivity contrast over space.


Introduction
In inverse scattering, we aim to retrieve the properties of objects from scattered waveforms. A common starting point is a nonlinear volume integral representation that relates the scattered field to contrasts in the medium properties with respect to a reference medium [1]. In acoustic inverse scattering, we can distinguish contrasts in compressibility and mass density. Although advanced algorithms have been developed for nonlinear inversion of both contrasts [2], it is still common practice to ignore the mass-density contrast and to rewrite the result as a nonlinear equation with a wave-speed contrast, which can be linearized under the Born approximation. Under this approximation, the wave-speed contrast can be retrieved directly from the data, for instance by filtered backpropagation as is routinely done in ultrasound diffraction tomography [3]. This methodology has found successful applications in ultrasonic breast imaging, given that both transmissions and reflections are available [4].
In reflection-mode diffraction tomography [5] or acoustic reflectivity imaging [6,7], the acoustic properties of an object are reconstructed from reflected waveforms only. Under these conditions, it appears difficult to retrieve the complete spatial spectrum of the contrast, especially when the object size is large compared to the wavelength and in the absence of low frequencies [8]. Since the acquired data are typically dominated by specular reflections, the retrieved image is often interpreted as the derivative of acoustic impedance along the object boundaries [9] (or interfaces). This seems somewhat inconsistent with the underlying theory, which is based on volume contrasts rather than on discontinuities at interfaces. In this paper, we propose a new formulation for reflectivity imaging that better suits the objective of quantifying the derivative of acoustic impedance along the object boundaries. In this formulation, we define reflectivity as the gradient of acoustic impedance in a medium with constant wave speed. We start our paper with a review of volume contrast imaging, which constitutes the core of reflection-mode diffraction tomography and reflectivity imaging. Then we derive a new form ulation for interface contrast imaging, which is designed to image acoustic interfaces with high resolution. Finally, we compare both imaging strategies by applying them to 2D synthetic ultrasound data of a cancerous human female breast.

Representation
Let x = (x 1 , x 2 , x 3 ) be a position vector in IR 3 and consider an unknown object that is irradiated from a variety of source locations at various frequencies. The object exists within a finite domain ID and is characterized by a spatially-dependent wave-speed contrast where c is the wave speed and c 0 is the (constant) wave speed outside ID. Throughout the medium, a constant mass density is assumed. We define p inc as the incident (pressure) field that would propagate through the homogeneous embedding and p as the total (pressure) field in the presence of the object. Let x R be a receiver location outside ID. The observed field at this receiver can be described by the Lippmann-Schwinger equation [1]: Here, γ = −jω/c 0 , where j is the imaginary unit and ω is the angular frequency. Further, G is the Green's function of the homogeneous medium with wave speed c 0 . This Green's function obeys the partial differential equation where ∇ denotes the gradient and δ(x) is a three-dimensional delta function. It is common practice to apply the Born approximation, p ≈ p inc , in the integrand, leading to This result can be used for modeling data at a specific frequency and receiver location, given the incident field at this frequency and knowledge of the volume contrast χ. Alternatively, the volume contrast can be retrieved from recorded data by linear inversion of equation (4).

Inversion
To utilize equation (4) for inversion, we assume that the domain can be represented by a family of piecewise-homogeneous subdomains (or voxels), which are indicated by index n. The integral over ID can be approximated by a finite summation of values of the integrand at x n , which we choose in the center of the subdomain. Assume further that data d = p − p inc are collected at various frequencies, source locations and receiver locations. The following data residual can be derived from equation (4): where χ n is the wave-speed contrast of subdomain n, ∆x is the mesh spacing, while j, r and s are indices for the frequency, receiver location and source location, respectively. We define the following error functional F (χ n ) = j,r,s ||r j,r,s (χ n ) || 2 j,r,s ||d j,r,s || 2 , which can be minimized with a conjugate gradient scheme [10]. This procedure is also known as Born inversion.

Interface contrast imaging
As we already mentioned in the introduction, Born inversion is relatively popular in ultrasound diffraction tomography, especially when both reflection and transmission measurements can be used. When the method is applied to reflection data only, it is often difficult to retrieve the full spatial spectrum of the volume contrast. Instead, the acquired image is typically interpreted as a measure for the derivative of acoustic impedance (i.e. the reflectivity). In this section, we propose an alternative formulation that better suits this interpretation. For the static field problem, the volume integral over a piecewise-homogeneous domain can be replaced by a collection of surface integrals enclosing the homogeneous subdomains, see e.g. [11]. This is possible because the static Green's function does not depend on the medium properties. If we assume a constant wave speed throughout the medium (i.e. c(x) = c 0 for x ∈ ID), the same strategy can be applied to the acoustic scattering problem.

Representation
To explain the methodology, we start with a simple scattering problem, where a constant acoustic impedance Z(x) = Z 1 is assumed throughout ID, while Z(x) = Z 0 in the exterior, see figure 1(a). The pressure field can be found by considering the boundary-integral formulation of this acoustic scattering problem, see section 8.2 of [1]. Further, we assume that the wave speed is constant throughout the medium and that the exterior of ID is strictly homogeneous with acoustic impedance Z 0 . For the pressure field at an arbitrary receiver outside ID, we have In the integrand of this expression, we find the particle velocity v and pressure p at the boundary, radiating towards the receiver via the Green's function G and its gradient ∇G, respectively. We proceed by eliminating the particle velocity from this result by subtracting another identity. To derive this identity, we apply the acoustic reciprocity theorem [1] to the inner domain ID with impedance Z 1 , while keeping the observation point x R outside ID. This procedure yields When we multiply this identity by Z 0 /Z 1 and subtract the result from equation (7), while using the continuity of the pressure and the normal component of the particle velocity, we find This equation can be linearized by assuming that p ≈ p inc in the integrand. However, we can improve this approximation as follows. We let the point of observation x R approach to an arbitrary point x on ∂ID and we approximate the integral in equation (9) by the contribution from the singular point x R = x. This leads to where we have defined the reflectivity factor Equation (11) has a similar form as the Kirchhoff approximation that is sometimes applied in seismic imaging [12]. We derived this approximation in a medium with constant wave speed rather than in a medium with constant mass density. The approximation can be used to express the (unknown) field p at ∂ID in terms of the (known) incident field p inc . Substitution into equation (9) yields Akin to the discrete representation that we discussed for volume imaging (i.e. equation (5)), we assume that the acoustic medium can be represented by a finite number of piecewisehomogeneous subdomains. More precisely, we replace domain ID by a family of disjoint subdomains ID n with boundaries ∂ID n , which are locally homogeneous and mutually separated by an infinitesimal distance h, see figure 1(b). We will now derive a relation akin to equation (13) for this configuration. Let ID be the union of all subdomains, for which equation (7) is valid and can be written as For each subdomain, we find an identity akin to equation (8): For each n, we multiply equation (15) with Z 0 /Z n and subtract the result from equation (14), yielding We define S as a collection of interfaces, which are placed in between adjacent subdomains (at a distance h/2), see figure 1(b). This collection also includes an outer boundary which encompasses the total domain ID (also at a distance h/2). As h → 0, the boundaries of adjacent subdomains align and their contribution can be merged. Therefore, the sum of integrals over boundaries ∂ID n can be replaced by a single integral over the interfaces in S. Consequently, equation (16) can be rewritten as where Z + (x) = Z(x + hν(x)) and Z − (x) = Z(x − hν(x)) denote the acoustic impedances at both sides of the interface at x ∈ S. In this result, we have used where ∇ R G is the gradient of the Green's function at x R . Since the right-hand side of equation (17) is invariant for the sign of ν , the direction of the normal of S is arbitrary and will be defined later. To obtain the field that is required in the integrand, we let the point of observation x R in equation (16) approach an arbitrary point x ∈ S and we approximate those integrals where this point is singular, yielding where we have defined a spatially-dependent directional reflectivity factor as When we substitute this approximation into equation (17), we find This result can be used to model data in a medium with constant wave speed, given the interface contrast R. Alternatively, the interface contrast can be retrieved by linear inversion of equation (21).

Inversion
In order to discretize our result, we assume a regular mesh of subdomains (indicated by index n) at x n (which we chose in the center of the subdomain). In each direction k ∈ {1, 2, 3}, we define an interface with midpoint x (k) n = x n + 1 2 ∆x i k , where i k is a unit vector of the coordinate system and ∆x is the mesh spacing. Let ν(x (k) n ) = i k be the normal of the interface, while S is the collection of all interfaces in all directions. The directional reflectivity of subdomain n in direction k can be expressed by discretizing equation (20) as To retrieve the reflectivity from the recorded data d j,r,s , we derive the following data residual after discretizing equation (21): We define an error functional as which can be minimized by a conjugate gradient scheme. In this way, the reflectivity contrast of each interface (specified by direction k) of each subdomain n can be retrieved.

Visualization
To visualize the reflectivity, the absolute value of each sample |R k,n | can be plotted as an interface with area (∆x) 2 and normal i k at its position x (k) n . By choosing an appropriate grayscale, we can distinguish between strong and weak reflectivity values. We can also visualize the polarity of the interface contrast with respect to some reference location x 0 (for instance the origin). For this purpose, we use the following definition to adjust the polarity of the reflectivity values: where sign k (x) denotes the sign of the kth dimension of x. By plotting the values of R 0 k,n with a polarity-dependent color, we can visualize whether a contrast is positive or negative with respect to the reference location x 0 .
The acoustic impedance can be estimated throughout the domain by integrating the retrieved reflectivity contrast over space. For this purpose, we write equation (22) as Since Z = Z 0 at the outer boundaries of ID is known, equation (26) can be applied recursively in an arbitrary direction i k to retrieve Z in the domain's interior.

Numerical example
To illustrate the theory presented in this paper, we apply numerical modeling and inversion to a 2D synthetic model that has been used in previous studies [13]. The model is derived from an MRI scan of a cancerous human female breast [14] and has been discretized into  on a refined grid with ∆x/5 spatial sampling using the free software package fdelmod [15]. For each source/receiver pair, we use 512 samples of data in the time domain with time sampling ∆t = 10 −6 s.

Volume contrast imaging
We start by evaluating volume contrast imaging. For this purpose, the error functional in equation (6) is minimized with a conjugate gradient scheme in order to obtain the wave-speed contrast throughout the computational domain. By allowing only real-valued updates, we enforce that the retrieved contrast is also real-valued. In figures 3(a) and (b), we show the retrieved contrast and its absolute value after one iteration (for visual purposes the inversion results have been interpolated to a grid of 2048 × 2048 pixels with ∆x/16 spatial sampling). Although the main features of the breast can be recognized (including the tumor), the resolution is relatively poor. Some improvements are observed after ten iterations, see figures 3(c) and (d). A higher number of iterations does not necessarily improve these results. This confirms the general consensus that the full spatial spectrum of the volume contrast can not be retrieved by Born inversion, given that only reflected waveforms are collected [8]. This can be attributed to the relative large size of the objects compared to the wavelength and the finite temporal frequency content of the data.

Interface contrast imaging
We use the same data and mesh to evaluate the performance of interface contrast imaging. For reference, we compute the directional reflectivity contrast at all pixel interfaces from the acoustic impedance model with the help of equation (22), see figures 4(a) and (b). We also compute R 0 with respect to the origin and the absolute value of the contrast, see figures 4(c) and (d).
Next, we retrieve these contrasts from the modeled data by minimization of the error functional in equation (24). In figures 5 and 6, we show the results after one and ten iterations, respectively. Note the superior resolution of these images compared to those of volume contrast imaging. The acoustic impedance can be estimated from the retrieved reflectivity contrast by applying equation (26) recursively. In figures 7(a) and (b), this is done after one iteration, using recursion in the x 1 -and x 2 -direction, respectively. This result improves significantly when the procedure is applied to the retrieved reflectivity values after ten iterations, see figures 7(c) and (d). Note that the acoustic-impedance contrast has been retrieved relatively well compared to the wave-speed contrast that was retrieved by volume contrast imaging.
To investigate the accuracy of these results quantitatively, we plot the normalized error norm E = log 10 (( j,r,s |r j,r,s | 2 )/( j,r,s |d j,r,s | 2 )) (on a logarithmic scale) for both methodologies as a function of iteration number, see figure 8. Note that interface contrast imaging converges faster than volume contrast imaging and leaves a smaller error norm at any number of iterations.

Discussion
Several alternative representations have been derived for imaging singular functions at reflecting interfaces, especially in the context of seismic imaging [16][17][18]. The general starting point for these representations is to consider constant mass density and a contrast in wave speed. In this paper, we reformulated the problem by considering constant wave speed and a contrast in acoustic impedance. As shown, this approach leads to a novel representation, which can be used for straightforward linear modeling and inversion in low-contrast media. We have tested this approach on a synthetic breast model with relatively large objects (compared to the wavelength) and a closed acquisition boundary. To avoid an inverse crime, the data have been computed by finite-difference modeling. The sensitivity to other types of noise is beyond our present scope. Also, the performance for different acquisition configurations and for more realistic biological tissues (with variations at a broader range of spatial scales) has not been considered. It is observed that interface contrasts are typically sparser than volume contrasts. In the numerical example that was discussed in this paper, for instance, only 3% of the interface contrast coefficients are non-zero, versus 23% of the volume contrast coefficients. Hence, the proposed representation in terms of interface contrast may also be used for compressive sensing, where one typically inverts for a sparse representation that fits the observed data by minimizing the number of non-zero coefficients [19].

Conclusion
In volume contrast imaging, we aim to retrieve a volume contrast in wave velocity under the assumption of constant mass density. The methodology is derived by linearizing the Lippmann-Schwinger equation under the Born approximation, leading to a linear relationship between data and contrast that can be solved by iterative inversion. Although this type of imaging has a good track record when both transmissions and reflections are acquired, it is difficult to retrieve the full spatial spectrum of the volume contrast from reflected waveforms only. When applied in pure reflection mode, we typically obtain an approximate image of object boundaries with relatively poor resolution.
In this paper, we propose a novel method to image acoustic interfaces with high spatial resolution. To achieve this, we have assumed that the medium can be partitioned into local subdomains. The acoustic impedance is assumed constant within each subdomain, while the wave speed is constant throughout all space. Under these conditions, a linear relationship can be established between the recorded data and a directional reflectivity contrast, which is nonzero only at a collection of interfaces between the subdomains. In interface contrast imaging, the reflectivity contrast is retrieved by iterative inversion. If necessary, the acoustic impedance can be estimated by integrating the retrieved reflectivity contrast over space.
We have applied volume and interface contrast imaging to synthetic pulse-echo data. For this study, we have used a 2D model of a human female breast, which contains realistic variations in both wave speed and acoustic impedance. With volume contrast imaging, we retrieve an image of the acoustic interfaces with relatively poor resolution. Interface contrast imaging enables us to improve the resolution considerably. Moreover, the methodology converges faster to a lower error norm, making it especially attractive for real-time imaging applications.