Nonlinear simultaneous reconstruction of inhomogeneous compressibility and mass density distributions in unidirectional pulse-echo ultrasound imaging

In diagnostic ultrasound imaging, the image reconstruction quality is crucial for reliable diagnosis. Applying reconstruction algorithms based on the acoustic wave equation, the obtained image quality depends significantly on the physical material parameters accounted for in the equation. In this contribution, we extend a proposed iterative nonlinear one-parameter compressibility reconstruction algorithm by the additional reconstruction of the object’s inhomogeneous mass density distribution. The improved iterative algorithm is able to reconstruct inhomogeneous maps of the object’s compressibility and mass density simultaneously using only one conventional linear transducer array at a fixed location for wave transmission and detection. The derived approach is based on an acoustic wave equation including spatial compressibility and mass density variations, and utilizes the Kaczmarz method for iterative material parameter reconstruction. We validate our algorithm numerically for an unidirectional pulse-echo breast imaging application, and thus generate simulated measurements acquired from a numerical breast phantom with realistic compressibility and mass density values. Applying these measurements, we demonstrate with two reconstruction experiments the necessity to calculate the mass density in case of tissues with significant mass density inhomogeneities. When reconstructing spatial mass density variations, artefacts in the breast’s compressibility image are reduced resulting in improved spatial resolution. Furthermore, the compressibility relative error magnitude within a diagnostically significant region of interest (ROI) decreases from 3.04% to 2.62%. Moreover, a second image showing the breast’s inhomogeneous mass density distribution is given to provide additional diagnostic information. In the compressibility image, a spatial resolution moderately higher than the classical half-wavelength limit is observed.

diagnostic information. In the compressibility image, a spatial resolution moderately higher than the classical half-wavelength limit is observed.
(Some figures may appear in colour only in the online journal)

Introduction
Ultrasound imaging modalities are applied as non-invasive, non-ionizing, real-time and costeffective diagnostic imaging techniques in several medical disciplines. As an alternative to x-ray mammography and magnetic resonance imaging (MRI), ultrasound examination of the female breast has been used in mammography diagnostics. Furthermore, due to the breast's direct spatial accessibility, automated ultrasound breast scanners and adapted algorithms, able to reconstruct different physical material parameters (e.g. speed of sound, attenuation), have been developed in the last few years (Duric et al 2005, André et al 2008, Hansen et al 2008. These technically sophisticated scanners mostly consist of two transducer arrays lying opposite in the same plane, or make use of a strong acoustic reflector across from a single transducer array. Data are acquired by rotating the transducer configuration around the breast, and reflected and/or transmitted ultrasound waves are detected for every fixed rotation angle and slice. Based on these measurements, advanced reconstruction concepts are utilized to image the breast's inner anatomy. The problem of reconstructing an object's inner structure from detected scattered measurements is called the inverse scattering problem. However, established ultrasound imaging techniques are solely able to reconstruct qualitative maps of variations in the tissue's acoustic impedance. Alternative concepts, applying a simplified linearized physical wave propagation model obtained by the Born or Rytov approximation, exist to reconstruct quantitative maps of the object's material parameters (Devaney 1985, Natterer andWübbeling 2001). Hence, an approximation of the scattered acoustic field is used to calculate a functional relationship between the Fourier transforms of the object's material parameters and the measurements, respectively. Applying this relationship, the object's material parameters are then reconstructed by the Fourier inversion. However, the approximation of the scattered acoustic field applying the Born approximation is just valid for weakly scattering objects and takes only single scattered waves into account. Due to these significant limitations in the wave propagation model, reconstructed images of objects with strong variations in the acoustic impedance providing multiple scattering (e.g. female breast) can be affected by artefacts (Simonetti et al 2008).
In order to avoid reconstruction artefacts obtained by ignoring multiple scattering in the wave propagation model, nonlinear diffraction tomographic reconstruction approaches have been published in several contributions (van den Berg 2002, Kowar 2005, van Dongen and Wright 2007, Natterer 2008a, Natterer 2008b, Natterer 2010a. In his contributions, Natterer proposes an iterative nonlinear algorithm for several scan configurations (one fixed linear transducer array, two fixed linear transducer arrays positioned opposing each other, and one fixed ring transducer array) able to reconstruct inhomogeneous maps of the tissue's compressibility under multiple scattering. The algorithm is based on the acoustic wave equation for inhomogeneous media having a space-varying speed of sound; these variations are solely assumed to be caused by a varying compressibility, whereas the mass density is considered to be constant (Morse andIngard 1986, Blackledge 2005). The reconstruction of compressibility is realized by solving the corresponding inverse scattering problem applying the iterative Kaczmarz method. Therefore, a numerical simulation of a forward wave propagation and a resulting residual backpropagation, using an explicit finite-difference timedomain (FDTD) scheme, is utilized to update the object's inhomogeneous compressibility distribution iteratively. The most simplified scan process is discussed in (Natterer 2010a) employing only one fixed conventional linear transducer array for wave transmission and detection. For this scan configuration, it is deduced from the Born approximation that low spatial frequency components of the object to be imaged can only be reconstructed if the excitation pulse contains low temporal frequencies (Natterer 2010a), or a strong reflector (e.g. metal plate) serving as an acoustic mirror is used (Natterer 2010b). Natterer transfers this result to the nonlinear compressibility reconstruction problem, and thus applies low temporal frequencies in the excitation pulse.
The assumption of a constant mass density distribution as considered by Natterer is widely-used when developing ultrasound reconstruction algorithms. However, this assumption certainly models a simplification of the actual tissue characteristics. This simplification could cause image artefacts if a reconstruction algorithm ignores mass density inhomogeneities while using measurements acquired from a tissue with significant spatial mass density variations. In order to avoid these artefacts, several algorithms, able to reconstruct both parameters simultaneously, have been proposed in various contributions: Devaney derives a generalized projection-slice theorem and filtered backprojection algorithms (Devaney 1985), and in (Anastasio et al 2005), an algebraic method in the Fourier space is proposed. However, both methods take only single scattered waves into account. Kowar uses an integral equation method to estimate the object's mass density and speed of sound, and discusses the reconstruction under multiple scattering for the case of a constant mass density and a space-varying speed of sound (Kowar 2005). In (van Dongen and Wright 2007), a contrast source inversion method for reconstructing both material parameters under multiple scattering is proposed. However, this method needs both the detected acoustic pressure and the particle velocity for reconstruction.
In this contribution, we develop a novel algorithm able to reconstruct inhomogeneous maps of the tissue's compressibility and mass density simultaneously under multiple scattering. Unlike to previous contributions, the algorithm uses only the detected acoustic pressure for reconstruction, and is based on an acoustic wave equation including spatial compressibility and mass density variations as well as the Kaczmarz method. We present in section 2 a mathematical derivation along the lines of the approach described in (Natterer 2010a) for this iterative nonlinear two-parameter reconstruction algorithm. In section 3, we investigate the potential performance of the improved algorithm numerically for an unidirectional pulse-echo breast imaging application, and evaluate the resulting image reconstruction quality in contrast to the quality achieved by the previously described compressibility reconstruction method.

Scan configuration and wave equation for inhomogeneous media
A schematic of the utilized two-dimensional unidirectional pulse-echo scan configuration is shown in figure 1. Denoting x 1 the lateral and x 3 the axial direction of a fixed linear transducer array, we consider an open simply connected medium = R 2 \ {x 3 0}. The transducer of width r 1 is located on the straight line {(x 1 , x 3 )| x 1 ∈ [−r 1 /2, r 1 /2], x 3 = 0}; therefore, the Sommerfeld radiation condition is fulfilled and the significant boundary part of is given by Let ⊂ be a closed inhomogeneous object to be reconstructed having a space-varying compressibility κ 1 and mass density ρ 1 , respectively; it is coupled to the transducer by the surrounding homogeneous medium \ (e.g. water) with constant compressibility κ 0 and constant mass density ρ 0 . For r = (x 1 , x 3 ) ∈ , the compressibility κ and the mass density ρ are thus defined by An incident sound wave propagating through the medium is then reflected at the medium's boundary as well as scattered at compressibility and mass density inhomogeneities within . The corresponding acoustic pressure p vanishes sufficiently fast at spatial infinity and satisfies the acoustic wave equation (Morse andIngard 1986, Blackledge 2005) as well as the initial condition where , ∇ and ∇· are the Laplacian, gradient and divergence operator with respect to r, andp denotes the second derivative with respect to t. The constant background speed of sound of both media is given by c 0 = (κ 0 ρ 0 ) −1/2 , and γ κ (r) = κ −1 0 κ(r) − 1 ∈ L 2 ( ) as well as γ ρ (r) = 1 − ρ 0 ρ(r) −1 ∈ L 2 ( ) represent relative compressibility and mass density deviations from the background values κ 0 and ρ 0 , respectively. Moreover, parameter T is the observation time. Note that the applied wave equation in (Natterer 2010a) is a special case of (2.1) using γ κ = f and γ ρ = 0, with f being the space-varying parameter to be reconstructed.
We consider the excitation of with an approximately cylindrical sound wave transmitted by a rectangular transducer element of width r 2 . This excitation is modelled by the boundary condition where ∂ n is the inner normal derivative, q describes the excitation pulse, a denotes the rectangular aperture function of width r 2 and s ∈ [(−r 1 + r 2 )/2, (r 1 − r 2 )/2] represents the centre position of the current transmitting source element. Summing up, the initial value problem (2.1)-(2.3) models a wave propagation through after medium excitation and is called the forward wave propagation.

Inverse scattering problem and Kaczmarz method
We consider a natural number m of exciting source elements at centre positions s 1 , . . . , s m ∈ [(−r 1 + r 2 )/2, (r 1 − r 2 )/2]. The corresponding inverse scattering problem is to reconstruct the material parameters γ κ and γ ρ from a given measurement data set G = {g 1 , . . . , g m }, g i = g i (r, t ), (r, t ) ∈ × [0, T ], representing the detected acoustic pressure while using the source element at centre position s i for excitation. Denoting γ (r) = (γ κ (r), γ ρ (r)) T ∈ L 2 ( ) × L 2 ( ), we define for each g i the nonlinear residual operator R i : with p being the γ-dependent solution of (2.1)-(2.3) while using the same source element at centre position s i for medium excitation. The introduced set L 2 ( ) × L 2 ( ) is called the space of material parameters. Due to the acoustic pressure's nonlinear dependence on the material parameter inhomogeneities γ, the nonlinearity of R i thus corresponds to the argument γ; furthermore, the operator R i maps the inhomogeneities γ to the residual error on the boundary . Consequently, γ ∈ L 2 ( ) × L 2 ( ) is reconstructed by solving the nonlinear homogeneous system System (2.4) can be solved by utilizing the Kaczmarz method (Natterer 1995, Natterer andWübbeling 2001). The related material parameter update for the source at centre position s i is given by represents an a priori filter for the residual g i − p. This filtered residual is then backpropagated by the adjoint operator R * i into the space of material parameters. Thus, update scheme (2.5) is similar to the filtered backprojection method (Kak andSlaney 2001, Buzug 2008) used as a standard reconstruction technique in CT.
In order to simplify the update scheme, we ignore the a priori filter in (2.5) replacing it by the identity operator as suggested in (Natterer 1995). Consequently, the material parameter update for the source at centre position s i using a simplified version of the Kaczmarz method becomes (2.6) The effect of ignoring the a priori filter can be reduced by choosing a suitable relaxation parameter as discussed in section 3. In linear seismic imaging, the above introduced a priori filter (R i (γ )R i (γ ) * ) −1 corresponds to a matrix operator (RR * ) −1 , with R being a forward seismic migration operator and R * its adjoint. Also in linear imaging, direct inversion is often not feasible due to the size and the bad-conditioned structure of the Hessian matrix RR * . Instead, the Hessian is often replaced by a diagonal or identity matrix as suggested above for the nonlinear case. Furthermore, several methods are proposed to approximate the inverse Hessian in a more accurate way; comprehensive overviews can be found in (Plessix and Mulder 2004, Wang et al 2009, Ren et al 2011.

Adjoint operator and iterative reconstruction algorithm
Employing the material parameter update scheme (2.6) to develop an iterative nonlinear twoparameter reconstruction algorithm, we primarily determine the adjoint operator R * i . Hence, we prove the following

to the source at centre position s i is given by
with p being the solution of the initial value problem (2.1)-(2.3), and z denoting the solution of the final value problem that tends sufficiently fast to zero at spatial infinity.
Proof. See appendix.
The introduced final value problem (2.8)-(2.10) also models a wave propagation through the medium similar to the initial value problem (2.1)-(2.3), however, the residual error g i − p is utilized for medium excitation. Hence, we call (2.8)-(2.10) the backward residual propagation.
Update (2.6) using (2.7) can obviously be computed for each source element at centre position s i . This application is called a sweep, and after computing the first complete sweep applying all sources randomly, one has to calculate additional ones to improve the image reconstruction quality iteratively. Thus, after determining a suitable initial approximation with γ κ = γ κ n,i−1 and γ ρ = γ ρ n,i−1 . (iv) Material parameter update: Use (i) and (iii) to compute the material parameter update for each image point r ∈ by  (e)). The breast is coupled by water (f) to a schematic of a fixed linear transducer array used for wave transmission and detection. Each white element is randomly selected for medium excitation, and all grey and white elements receive the scattered echoes during the process of reconstruction. The black elements are inactive.
In order to improve convergence, we select different α -values α κ and α ρ in (2.11) and (2.12), respectively, and make them space-dependent as discussed in section 3. The introduced reconstruction algorithm can also be implemented for modified scan configurations (e.g. two fixed linear transducer arrays opposing each other, one fixed ring transducer array, etc). Moreover, the presented mathematical derivation and the resulting reconstruction algorithm are still valid for the three-dimensional case. Note that in case of γ ρ = 0, a sweep (a)-(d) is equivalent to the sweep derived in (Natterer 2010a), except for the neglected prefactor c −2 0 .

Numerical breast phantom
The self-created numerical breast phantom (300 × 500 grid points) together with a schematic of a fixed linear transducer array we used for our reconstruction experiments is illustrated in figure 2. In figure 2, the same numerical breast phantom as in (Hesse and Schmitz 2012) is shown; in contrast to this contribution, we reduced the number of transmitting elements from 100 to 20 and displayed several inactive elements to clarify the upper boundary condition. The phantom had a lifelike anatomy and size, and we utilized realistic compressibility and mass density values for varying tissue types and water as specified in table 1. We considered a transducer array with 100 active elements (grey and white elements) coupled to the breast by water; the rectangular elements had a width of r 2 = 1.5 mm ( = 4 grid points), and the kerf was 0.5 mm ( = 2 grid points). Moreover, we selected a limited number of m = 20 equidistant elements (white elements) for medium excitation. Each of these elements was used to transmit an approximately cylindrical sound wave, and all active elements were used to receive the scattered echoes. Note that for an experimental scanner construction, the patient would be lying facedown on an examination table with the breast immersed in a water tank and the transducer array fixed below the breast. In this case, the illustrated scan configuration in figure 2 needs to be rotated through 180 degrees. Table 1. Speed of sound, compressibility, mass density, relative compressibility deviation and relative mass density deviation values for varying tissue types and water (Foster et al 1984, Duck 1990

Numerical implementation
In order to execute steps (i) and (iii) of the above introduced iterative sweep (a)-(d), we employed an explicit FDTD scheme as presented below. Hence, we approximated all derivatives in (2.1) and (2.8) by central differences with second-order error terms truncated; furthermore, we made use of a forward difference formula, again with second-order error terms truncated, to discretize (2.3) and (2.10) (Hirsch 2007). We applied the following time shifted modification of the proposed excitation pulse in (Natterer 2010a) with ω c = 2π f c being the angular centre frequency. The instant of time with maximum excitation energy supply is realized at t = t 0 and hence, it is essential to put the fixed value t 0 > 0 sufficiently large. The fractional bandwidth of the introduced excitation pulse is 83%. In (3.1), the excitation pulse amplitude A is normalized to 1. In order to obtain typical physical pressure values, much higher amplitudes have to be used (e.g. A = 5 × 10 8 to achieve a peak negative pressure of about 0.2 MPa). It has to be taken into account that with an unnormalized amplitude A, the space-dependent relaxation parameters as presented below have to be divided by A 2 according to the linearity of the FDTD scheme.
Applying an analogous calculation to that in (Natterer 2008b), one can show for the linearized wave propagation model achieved by the Born approximation that low spatial frequency components of γ κ and γ ρ can only be reconstructed if the excitation pulse q contains low temporal frequencies. We transferred this result to our extended reconstruction approach, and thus initially chose f c = 0.1 MHz with t 0 = 17 μs to preestimate γ κ and γ ρ by calculating some sweeps (a)-(d). We subsequently employed the resulting material parameter update as a new initial approximation to improve the reconstruction quality of γ κ and γ ρ with (a)-(d) using f c = 0.2 MHz with t 0 = 9 μs. Due to that choice of frequencies and the Courant-Friedrichs-Lewy condition (Inan and Marshall 2011), we put x 1 = x 3 = 0.5 mm and t = 0.2 μs in the corresponding FDTD scheme using it for computations with both centre frequencies. The FDTD scheme and the medium's excitation related to the forward wave propagation are given by withc = c 0 t/ x 1 being the Courant number, i denoting the discrete time step as well as j and k representing the discrete spatial steps in x 1 -and x 3 -direction, respectively. The backward residual propagation scheme calculating z is identical to (3.2), whereas (3.3) has to be modified using (2.10) to implement the residual error excitation. The initial value problem (2.1)-(2.3) in sweep step (i) can be solved directly with (3.2) and (3.3), whereas a time reversal (Natterer and Wübbeling 2001) is needed to solve (2.8)-(2.10) in step (iii). In order to model absorbing boundaries, we implemented second-order Mur boundary conditions in all except for the first, second, penultimate and last grid point on the phantom's left, right and lower boundary; furthermore, we used first-order Mur boundary conditions for the second and penultimate grid point (Inan and Marshall 2011). Moreover, we applied the composite trapezium rule to calculate approximations of the correlation integrals in (2.11) and (2.12).

Reconstruction setup
In order to validate our derived reconstruction algorithm, two numerical reconstruction experiments were conducted. Initially, we acquired two synthetic measurement data sets G 1 = {g 1 1 , . . . , g 1 20 } and G 2 = {g 2 1 , . . . , g 2 20 } as previously introduced, using (3.2) and (3.3) with f c = 0.1 MHz and f c = 0.2 MHz, γ κ and γ ρ as listed in table 1 and T = 179.8 μs. Data set G 1 was primarily employed to reconstruct low spatial frequency components (spatially extended homogeneous regions) of the utilized numerical breast phantom, whereas G 2 was used for higher spatial components (sharp tissue boundaries). Both numerical experiments applied the same data sets G 1 and G 2 for reconstruction.
In the first reconstruction experiment, we evaluated the induced image artefacts caused by ignoring spatial mass density variations during the process of reconstruction. Hence, we considered the compressibility to be the only space-varying material parameter, and thus executed each sweep (a)-(d) with γ ρ = 0. More precisely, we chose γ κ = 0 as initial approximation, and applied G 1 to preestimate γ κ . We utilized the resulting preestimation as a new initial approximation to reconstruct γ κ finally using G 2 . For both estimations, we chose α κ (r) = 10 10 x 3 m −1 proportional to the perpendicular transducer distance as suggested in (Natterer 2010a), and computed 1000 iterative pre-and mainsweeps. The executed reconstruction was equivalent to that in (Natterer 2010a), however, the applied measurement data sets G 1 and G 2 took both compressibility and mass density inhomogeneities into account.
In the second reconstruction experiment, we considered both material parameters to be space-varying and did their reconstruction simultaneously. We applied γ κ = γ ρ = 0 as initial approximation to preestimate γ κ and γ ρ using G 1 , and employed the final results to improve the image reconstruction quality utilizing G 2 . Again, we calculated 1000 iterative pre-and mainsweeps for both estimations, however, we put α κ (r) = 10 11 x 3 m −1 and α ρ (r) = 10 10 x 3 m −1 .

Reconstruction results
The compressibility reconstruction image and the corresponding relative error magnitude achieved in the first experiment are shown in figure 3. Water, connective, fat, muscle and tumourous tissue are correctly placed (a). The reconstruction quality is quantitatively satisfying and the mean value within the indicated diagnostically important ROI amounts to 3.04% (b). However, as can be seen from (a), the spatial resolution in the breast's diagnostically significant inner part is moderate. Glandular tissue is not sharply bounded from the surrounding connective tissue and thus, the reconstruction's diagnostic value is reduced. Indeed, the tumour of about . Compressibility reconstruction image and related relative error magnitude obtained in the first experiment ignoring spatial mass density variations during reconstruction. Compressibility reconstruction after 1000 pre-and mainsweeps ranging from 3 × 10 −10 to 6.5 × 10 −10 Pa −1 (a), and related relative error magnitude ranging to 15% (b). Compressibility values on the red line are used for spatial resolution comparison (a). Relative error magnitude values within the blue diagnostically important ROI are applied to determine the quantitative reconstruction quality (b). 4.5 mm diameter is observable, but, as expected from the introduced compressibility values for connective and tumourous tissue in table 1, it has a low contrast to the surrounding connective tissue.
The simultaneous reconstruction results of both material parameters compressibility and mass density as well as related relative error magnitudes obtained in the second experiment are displayed in figure 4. In each reconstruction image, water and all tissue types are correctly located (a), (c). As is apparent from (b) and (d), both reconstructions are quantitatively superior in contrast to the first experiment; the mean values in the diagnostically important ROIs amount to 2.62% (b) and 1.55% (d), respectively. In contrast to the first experiment, the spatial resolution in the breast's inner part is improved, and glandular tissue is sharply bounded from the surrounding connective tissue in both images (a) and (c). Again, the tumour has the expected low-contrast resolution to connective tissue in the compressibility image (a). However, in the mass density image (c), the tumour is clearly visible having a high contrast to connective tissue; this corresponds to the introduced mass density values for connective and tumourous tissue in table 1.
Comparisons of spatial material parameter resolutions based on compressibility and mass density values calculated in the first and second experiment as well as values from the original numerical breast phantom are shown in figure 5. As evident from the left and right side of (a), both algorithms reconstruct compressibility values of water and fat tissue almost correct. In order to reconstruct connective and glandular tissue in the breast's diagnostically important inner part, our extended algorithm applied in the second experiment seems to be able to detect minor compressibility variations using more oscillations for approximation in contrast to the algorithm used in the first experiment. Accordingly, our algorithm utilized in the second experiment is thus able to reconstruct sharp boundaries between the connective and the glandular tissue resulting in improved spatial resolution. As stated in (b), our algorithm is additionally able to reconstruct the mass density almost correctly. Water, connective, fat and glandular tissue values in the left and middle part are quantitatively satisfying, however, the approximation of connective and fat tissue values on the right side is less good.
To investigate the spatial resolution of reconstruction results in figure 5 quantitatively, we determined the transfer function of the reconstruction method. For this, we calculated the one-dimensional Fourier transforms of the input material parameter distribution and the  . Compressibility and mass density reconstruction images and related relative error magnitudes obtained in the second experiment providing spatial mass density variations during reconstruction. Compressibility reconstruction after 1000 pre-and mainsweeps ranging from 3 × 10 −10 to 6.5 × 10 −10 Pa −1 (a), and related relative error magnitude ranging to 15% (b). Mass density reconstruction after 1000 pre-and mainsweeps ranging from 900 to 1100 kg m −3 (c), and related relative error magnitude ranging to 15% (d). Compressibility and mass density values on the green lines are used for spatial resolution comparison (a), (c). Relative error magnitude values within the blue diagnostically important ROIs are applied to determine the quantitative reconstruction qualities (b), (d). reconstructed material parameters' distribution. The −6 dB spatial cut-off frequencies giving the upper resolution limit of the reconstruction algorithms were calculated. For κ, the cut-off frequency was 0.42 mm −1 in the first experiment and 0.49 mm −1 in the second experiment, demonstrating the improved resolution. For ρ in the second experiment, a cut-off frequency of 0.24 mm −1 was calculated. The theoretical half-wavelength limit (Wolf 1969) for a frequency of 0.283 MHz (upper band-limit of the excitation pulse) and the slowest speed of sound of 1437 m s −1 is 0.39 mm −1 . For κ, this limit is exceeded. This effect may be caused by multiple scattering as discussed in (Schatzberg and Devaney 1992, Simonetti 2006and Simonetti et al 2008. Simonetti demonstrates the possibility of superresolution by evanescent waves generated by scatterers that are scattered again, encoding subwavelength information into the scattered far-field. Summing up, in order to avoid reconstruction artefacts in the compressibility image in case of tissues with significant mass density inhomogeneities, it is essential to take spatial mass density variations into account during the process of reconstruction. As is apparent from the second reconstruction experiment, the mass density image potentially provides the detection of abnormal tissues (e.g. tumourous tissue) having a low-contrast compressibility resolution, but a high-contrast resolution in the mass density image.

Discussion and future work
The above presented numerical results demonstrate the potential reconstruction performance of our extended algorithm for an unidirectional pulse-echo breast imaging application. We apply a conventional linear transducer array at a fixed location for data acquisition, instead of rotating a complex transducer configuration around the object to be imaged. Unlike to previous studies, the derived iterative nonlinear approach is able to reconstruct the object's inhomogeneous compressibility and mass density separately under multiple scattering using only the detected acoustic pressure. Several published algorithms are solely able to reconstruct inhomogeneous maps of the object's speed of sound as a composition of the unknown fundamental material parameters compressibility and mass density, reconstruct inhomogeneous compressibility distributions under the assumption of a constant mass density, or apply both the detected acoustic pressure and the particle velocity for calculation. However, the advantage of simple data acquisition is reduced by the requirement to use a low-frequency excitation pulse. Moreover, this excitation pulse as well as the constant background speed of sound of both media have to be known precisely for reconstruction.
Obviously, due to the algorithm's complexity and the large number of calculated sweeps, reconstruction at video rate is currently not feasible. Indeed, acceptable image reconstruction results can be obtained by using a clearly reduced number of sweeps. However, high-resolution benefits and the diagnostic value are reduced in these images. In order to accelerate the image reconstruction process, the explicit FDTD scheme of both the forward wave propagation and the residual backpropagation can be effectively parallelized on parallel computing platforms (e.g. graphic processing unit (GPU)).
In future studies, the derived algorithm needs to be implemented on an appropriate parallel computing platform to accelerate the image reconstruction performance. Moreover, the approach has to be evaluated for tissues with attenuation. Most likely, the reconstruction quality will be reduced for deep-seated image points in case of significant attenuating objects. In this case, a suitable extension of the acoustic wave equation including attenuation has to be applied to improve the introduced algorithm by the additional reconstruction of an inhomogeneous attenuation map. Another component that needs to be investigated is the impact of modified space-dependent relaxation parameters on the image reconstruction quality and speed of convergence. Furthermore, the potential superresolution properties of the reconstruction method have to be investigated systematically, and the reconstruction performance of the derived approach needs to be verified experimentally.

Conclusion
In the present contribution, we extended a proposed iterative nonlinear one-parameter compressibility reconstruction algorithm by the additional reconstruction of the object's inhomogeneous mass density distribution. The improved reconstruction approach is able to reconstruct inhomogeneous maps of the tissue's compressibility and mass density simultaneously applying the Kaczmarz method for iterative material parameter update. We validated the performance of our algorithm numerically, and demonstrated for an unidirectional pulse-echo breast imaging application using only one fixed conventional linear transducer array for wave transmission and detection, the necessity to provide spatial mass density variations during reconstruction for tissues with significant mass density inhomogeneities. Reconstructing spatial mass density variations, artefacts are qualitatively and quantitatively reduced in the compressibility image resulting in improved spatial resolution. Moreover, a second diagnostic image, showing the inhomogeneous map of the tissue's mass density distribution, is additionally given to exam tissues from an alternative diagnostic point of view. Furthermore, the achieved spatial resolution of compressibility reconstruction images was observed to moderately exceed the classical half-wavelength resolution limit, an effect that has to be investigated thoroughly by simulations and experiments. According to the obtained numerical results, our extended approach has the potential to reconstruct quantitative highresolution compressibility and mass density images and to outperform the image reconstruction quality achieved by established ultrasound imaging techniques.
Although not discussed in detail in this contribution, the mathematical derivation and the resulting reconstruction algorithm can easily be generalized to three dimensions. Hence, three-dimensional inhomogeneous compressibility and mass density distributions can be reconstructed using a fixed two-dimensional planar transducer array for wave transmission and detection. Furthermore, the algorithm can also be adapted for modified scan configurations (e.g. two transducer arrays, ring transducer, etc), even in a three-dimensional setup.

Appendix
In this appendix, the proof of theorem 1 is presented.