Phase-insensitive versus phase-sensitive ultrasound absorption tomography in the frequency domain

The sensitivity of phase-sensitive detectors, such as piezoelectric detectors, becomes increasingly directional as the detector element size increases. In contrast, pyroelectric sensors, which are phase-insensitive, retain their omni-directionality even for large element sizes, although they have significantly poorer temporal resolution. This study uses numerical models to examine whether phase-insensitive detectors can be used advantageously in ultrasound tomography, specifically absorption tomography, when the number of detectors is sparse. We present measurement models for phase-sensitive and phase-insensitive sensors and compare the quality of the absorption reconstructions between these sensor types based on relative error and image contrast metrics. We perform the inversion using synthetic data with a Jacobian-based linearized matrix inversion approach.


Introduction
Ultrasound tomography (UST) is an emerging approach for tomographic reconstructions which functions through the measurement of ultrasound waves after interaction with a target of interest.UST has found particular interest in the imaging of soft tissues, such as breasts, where ultrasound transmission tomography has been employed to reconstruct the absorption and sound speed profiles of the interior tissue [1,2,3,4], since cancerous breast tissue is known to have different absorption and sound speed properties compared to healthy breast tissue [5,1].An important aspect of assessing the performance of UST systems is the detector type, its size and its directional response.In particular, we consider here the possibility of using pyroelectic detectors which are phase-insensitve [6], versus the better known piezoelectric detectors which are phase-sensitve [7].We will here consider ultrasound transmission tomography for the purpose of acoustic absorption reconstructions only, as it best enables the comparison of PS and PI detector types in the parallel array geometry, which further allows for comparisons with traditional x-ray tomography approaches.PI sensors have previously been studied in the context of UST absorption reconstructions with the use of pyroelectric ultrasound sensors [6,8], which is the sensor type that we will use as a reference to model PI sensors in this paper.The key aspect of PI sensors which may prove beneficial in certain scenarios is their much flatter directional response curve due to the lack of phase-cancellation [6,9], which for PS sensors can lead to false absorption measurements.This is especially true in larger sensors, which have better signal-to-noise ratios but also suffer from strong directionality in PS sensors.
We evaluate the reconstruction qualities of PS and PI sensors for a parallel array geometry with a varying range of array rotation angles and number of source frequencies, for two choices of sensor width.The quality evaluation is done through two different contrast metrics: the modulation transfer function (MTF) [10,11], and the root-meansquare (RMS) contrast [12].The MTF provides a contrast measure as a function of spatial frequency in the reconstructions, whereas the RMS contrast provides a global contrast measure that is independent of spatial frequencies.
In section 2 we outline the inverse problem that we are solving in this paper, along with the source-sensor geometry being considered.Then, in section 3 we describe the models used for the forward problem, consisting of the acoustic field simulation and the detector measurement models.Section 4 covers the theoretical description of the image reconstruction, with a derivation of the Jacobian operator used in our linearised reconstrcution approach.Section 5 explains our numerical simulation setup, how the image quality is analysed, followed by the results for our comparison of PS and PI sensors with a focus on sensor size and data sparsity and their affect on reconstruction quality for the two sensor types.Finally, in section 6 we discuss the results and their implications for UST and future research, in particular the scenarios in which PI sensors may be able to produce better absorption reconstructions than PS sensors.

Inverse problem
Our goal is to reconstruct the acoustic absorption profile of a region located between an array of ultrasound sources and an array of sensors through transmission UST, using the parallel array geometry shown in Figure 1.This source-sensor geometry can be compared to previous work in both UST as well as traditional x-ray computed tomography.The parallel array consists of a linear array of N point sources S ωθ,n , n ∈ {1, . . ., N} opposing a linear array of M sensors taken to be simply the integral over a finite support function χ θ,m ⊂ χ, m ∈ {1, . . ., M} of width d.We consider an infinite Euclidean domain χ = R 2 , so there are no boundary conditions affecting the sound waves.The source and sensor arrays are rotated with respect to their centre point over a range of angles θ ∈ Θ, and the sources are driven at a range of frequencies ω ∈ Ω in order to capture the full data for the reconstruction problem.The size of the forward model matrix scales linearly with each of the number of sources, sensors, angles, and frequencies.
Figure 1: Parallel array source and sensor geometry where for a rotation of θ the m-th sensor occupies a region χ θ,m ⊂ χ = R 2 .Each sensor in the array has a diameter of d and is located at a distance D from the source array.

Acoustic model
The forward model used in this paper is an acoustic Helmholtz equation with an absorption term implemented through a complex wavenumber, given by the family of equations which are parametrised by the source frequency ω, angle θ, and source index n.
Parameters which are controlled are written as indices, whereas the absorption and sound speed parameters, τ and c, which are defined by the medium are written in the function argument.Each constant frequency source S ωθ,n (x) gives rise to a complex acoustic pressure solution P ωθ,n (x; τ, c) through the differential operator L ω (x; τ, c) which has the form The absorption term τ is dimensionless, but can be related to the usual dimensional absorption coefficient α with units of dB cm −1 by the relation which is derived by defining the absorption coefficient to be the imaginary part of the complex wave number in (2).This form for the absorption term assumes a linear power law with respect to frequency, α = α 0 ω, where α 0 = τ /c.

Symbol
Jacobian operator of F at τ 0 , its discrete-continuous components, and discretised components f {i} ({x}; {a}) function f with variables {x}, experimental parameters {i}, and medium parameters {a} Table 1: Definitions for symbols used throughout this paper.
In this paper we are only interested in comparing absorption reconstructions, so we assume that the sound speed c is known and has a constant value throughout the medium χ for simplicity.The model does allow for heterogeneous sound speed and absorption.

Measurement model
In this paper we model the output of each PS sensor as proportional to the integral of the acoustic pressure over the sensor area, for each sensor in the array, and the output of each PI sensor as proportional to the integral of the squared acoustic pressure amplitude over the sensor area.These choices are designed to approximate the behaviour of piezoelectric (PS) sensors, which respond to pressure, and pyroelectric (PI) sensors, which respond to heating.The use of the squared pressure amplitude to model the pyroelectric sensor response follows from the observation that the directional response of a pyroelectric sensor correlates strongly with the directionality of the heat deposition in the sensor [9].Ultrasonic heat deposition is commonly taken to be proportional to the acoustic pressure amplitude squared [13].
In general we consider the measurement y ωθ,nm at the m-th sensor from the n-th source at frequency ω and array angle θ to be modelled as the composition of three operators as follows, where I θ,m is a sampling operator for the m-th sensor with rotation θ about the centre, M is a pointwise field transformation dependent on the type of sensor, and L −1 ω (x; τ, c) is the inverse Helmholtz operator for an absorption map τ .These operators act as follows: where 1 χ θ,m is the indicator function on the m-th sensor region at rotation θ.

Forward problem
We define a forward mapping where the absorption maps are described by square-integrable functions over χ, τ ∈ L 2 (χ), and the space of data is the complex space Y = C |Ω||Θ|N M .The components of the data vectors y ∈ Y are related to the absorption maps in L 2 (χ) by ( 4).Note that for PI sensors the imaginary part of the measurement is always zero.

Linear reconstruction scheme
Our reconstructions use a linear inversion scheme using the Fréchet derivative of the measurement model with respect to the parameter of interest, τ [14,15,16].The modelled measurement at a desired absorption distribution τ is related to the modelled measurement at a given absorption distribution τ 0 through the Taylor expansion of the modelled data at τ 0 : where h(x) = τ (x)−τ 0 (x).The first order linear operator the Fréchet derivative for our system evaluated at τ 0 , i.e. the linearization of the forward mapping from equation ( 8).In the continuous-discrete setting this can be called the Jacobian operator which we can define explicitly as follows: where i indexes over the full set of data parametrized by ω, θ, n, and m, and the column index is given by the continuum of positions x ∈ χ.We utilize the adjoint state method to compute the Jacobian operator in equation ( 10) one row at a time by taking the variational derivative of the modelled measurement y ωθ,nm (τ, c) with respect to the absorption parameter τ .We begin by taking an arbitrary variation in the absorption τ (x) → τ (x) + δτ (x) which results in a variation in the pressure P ωθ,n (x; τ ) → P ωθ,n (x; τ ) + δP ωθ,n (x; τ ), where δP ωθ,n (x; τ ) = P ωθ,n (x; τ + δτ ) − P ωθ,n (x; τ ) = ∂P ωθ,n (x;τ ) ∂τ δτ (x), and apply this variation to the modelled measurement: = (I θ,m • M)P ωθ,n + I θ,m M ′ (P ωθ,n )δP ωθ,n + O(δP 2 ωθ,n ) Ignoring higher order terms, we then have that the difference is given by δy ωθ,nm (τ, c) = y ωθ,nm (τ + δτ, c) − y ωθ,nm (τ, c) where we have used the relation The left argument of the inner product on the penultimate line of ( 12) can be understood as the adjoint field Z ω,nm through the adjoint equation where the right-hand side of ( 14) is the adjoint source term.The partial derivative of the Helmholtz operator L ω with respect to the absorption parameter τ is the imaging condition term, which can be evaluated to be Once the Jacobian has been computed at some absorption value of τ 0 (x) we can use standard matrix inversion schemes to solve for the absorption τ (x) through the linear approximation equation from Eq. ( 9) which defines the linearisation of our forward model from (8).For the inverse problem we have data vector g instead of the modelled data at the true absorption profile, y(τ ), so we define our inverse problem through the minimization of the residual in ( 16) along with a first order Tikhonov regulariser term added to better deal with the ill-posedness, ĥ = arg min where g is the measured data, y(τ 0 ) is the modelled data, ĥ is the optimal reconstructed difference in absorption maps τ − τ 0 , η ≥ 0 is the regularization parameter, and ∇ x , ∇ y are the spatial gradients in the x-and y-directions, respectively, in the domain χ.
If we discretize the domain χ to have L x pixels in the x-direction and L y pixels in the y-direction, so that we have a finite collection of position coordinates {x j } LxLy j=1 , then we can express the Jacobian operator from (10) as a (|Ω||Θ|NM) × (L x L y ) matrix with discrete components, ( J F | τ 0 ) ij .We have chosen to use MATLAB's built-in least-squares solver, LSQR, through the augmented matrix A and data vector b given by: to solve the minimization problem in Eq. ( 17).After reconstruction we further perform a frequency domain filtering step on the reconstruction ĥ(x) by removing very high spatial frequencies from the reconstruction which may arise due to small singular values in the augmented matrix A. We use the L-curve method to help choose the optimal regularization parameter η for each reconstruction.

Numerical experiments
Our modelled absorption profile is homogeneous with a constant dimensionless absorption of τ 0 (x) = 0.003, whereas the true absorption profile shown in figure 2 has an equal background value of τ (x) = 0.003 for x outside of the square, and a higher value of τ (x) = 0.006 for x within the square region.The sound speed is set to a constant c = 1540 m s −1 for both the modelled and measured data, which is a typical sound speed found in human tissue [17,18].Hence, we are assuming that the sound speed and background absorption are known exactly.From the absorption transformation equation (3) we can see that these dimensionless absorption values correspond to α(x)/2πω = 1.06 dB cm −1 MHz −1 for the background, and α(x)/2πω = 2.13 dB cm −1 MHz −1 for the square target, which are comparable to the absorption values measured in human tissue [17,19,20].
The simulation domain is a 40 mm × 40 mm square which is discretised into a 256 × 256 pixel grid, so we have that L x = L y = 256 and dx = dy = 0.15625 mm.The square target region in the true absorption profile τ has a side length of 50 pixels or 7.8125 mm, spanning the range [100, 150] × [150, 200] in pixels or [−4.4706, 3.3725] × [3.3725, 11.2157] mm in the spatial domain χ.The perfectly matched layer (PML) is implemented as a quadratic absorbing function [21] spanning 5 pixels (≈ 0.78 mm) on all sides of the domain.The low-pass filtering of the reconstructions is done using a circular cutoff in the frequency domain, where our chosen cutoff radius was set to remove spatial frequencies above 1.75 mm −1 , which corresponds to a circle with a 70 pixel radius in the 256 × 256 pixel frequency domain.
Our parallel array consists of 10 point sources spanning a distance of 30 mm and an array of 10 sensors of width d spanning the same 30 mm distance.The source and sensor arrays are 30 mm apart.The point sources each have an amplitude of 1 Pa, although this choice is arbitrary in the numerical setting as it only affects the scale of the measurements.
Our choices of τ 0 and c exhibit rotational symmetry, since they are both spatially constant, and hence we only need to compute the rows of our Jacobian using (13) for a single angle θ ∈ Θ.The rest of the rows can be computed by taking the sensitivity map from the computed row i ≡ ω, θ, n, m, and rotating it by an angle θ ′ −θ about the center of the parallel array to get the corresponding sensitivity map for row i ′ ≡ ω, θ ′ , n, m.This reduces the number of forward (P ) and adjoint (Z) field computations by a factor of |Θ|.For mediums with translation and reflection symmetry, in the case of identical sources and sensors, further optimisations may be made by computing a sensitivity map between source n and sensor m, and then performing the appropriate translation and reflection to get the corresponding sensitivity map between source n ′ and sensor m ′ that have the same, up to a reflection, relative displacement from each other as the original pair.This latter optimisation was not done for our computations, but it could further reduce the number of forward and adjoint field computations by a factor of N.
We want to compare the quality of the absorption reconstructions as three parameters are varied for each of the two sensor types: number of angles, number of source frequencies, and sensor size.The number of angles and frequencies control for the amount of data provided by the measurements, and hence give us information about how the reconstruction qualities scale with varying levels of data sparsity.The sensor size affects the directional response of the detector as well as the signal strength, and thus varying this parameter provides us with information about how well each sensor Table 2: Model parameters used for simulations.
responds to an increase in sensor size.
The sets of angles are constant increment subsets of θ ∈ [0 • , 180 • ), for the increments ∆θ = 7.5 • , 15 • , 30 • , and 60 • .The source frequencies are an odd number of constant increment frequencies centered around 2 MHz; we consider the two sets of frequencies Ω = {2} MHz and Ω = {1.5, 1.75, 2, 2.25, 2.5} MHz.Sensors of width 1 mm and 5 mm are considered, arranged into a linear array of 10 sensors spanning a distance of 30 mm.The 5 mm sensor array contains overlapping sensors, which can be realised physically by spatially translating a non-overlapping sensor array.Key simulation parameters are listed in table 2.
Both noise-free and noisy data are used for reconstructions.For noisy data, additive Gaussian noise scaled by a factor of 1% of the maximum sensor signal has been applied.

Contrast analysis
In order to quantitatively judge the quality of the reconstructions, we use two different contrast metrics: the modulation transfer function (MTF), and the weighted root mean square (RMS) contrast weighted by the maximum value of the reconstruction.See Appendix A for more details on these contrast metrics.To make the computation of these contrast metrics as easy as possible, the reconstruction target under consideration here is a square of constant absorption, shown in Figure 2.
The MTF FWHM for the true absorption difference profile is infinite, since the edge is a perfect step function in theory.The maximum weighted RMS contrast is 0.5 for the true absorption difference profile in a region of interest where half the region consists of the background profile, and the second half consists of the target profile, as indicated by the dashed-line box in figure 2.

Results
We find that in the noise-free cases the PS sensors are typically dominant in both MTF FWHM and RMS contrast measures except for some of the 3-angle reconstructions.When 1% noise is added, the PI sensors perform closer to the PS sensors in these contrast metrics, and sometimes outperform them.We can see examples of these sparse data reconstructions for 5 mm sensors in figure 3.In the full data regime, where we have 24 angles and 5 frequencies, the PS sensors beat the PI sensors for the MTF FWHM contrast metric for both noiseless and noisy data cases, but PI sensors can still achieve superior RMS contrast when using the smaller 1 mm sensors.The full data reconstructions for the 1 mm sensors can be seen in figure 4.
The contrast metrics for the different cases are summarised in figures 5 and 6.The greatest advantage for PI sensors is achieved with noisy data when the sensor is large in the sparse data regime, which we can see visually in figure 6c and 6d.

Conclusion
We have shown that phase-insensitive ultrasound sensors can outperform traditional phase-sensitive sensors in limited data situations, especially with large sensorsproducing superior contrast in absorption reconstructions both globally and across different spatial frequencies.Some uncertainty does arise from the choice of regulariser as well as regularisation parameter for the reconstructions.
We chose a standard Tikhonov regulariser that penalises higher magnitude gradients, which does a good job of making the reconstructions stable especially with the presence of noise.However, the amount of  regularisation needs to be carefully balanced in order to make fair comparisons across different reconstructions.Greater regularisation smooths out the reconstruction which improves global contrast at the cost of losing edge sharpness, and hence reducing the MTF FWHM.Similarly, lower regularisation leaves more artifacts in the reconstruction which decrease contrast, but generally leaves the edge profile more sharp, leading to a better MTF FWHM measure.The L-curve method was quite good at picking a regularisation parameter that balanced these factors quite well, but fine tuning through visual inspection was necessary in many cases to find the ideal parameter.
There is also a possibility that our simple linear inversion approach leaves one of the sensor types at a comparative disadvantage in certain situations.Nonlinear inversion schemes could be tested in select cases to see if notable relative differences in these contrast measures arise.In each case there are two line pairs distinguished by the markers along the curves, corresponding to the noise-free (l) and noisy (♦) cases.

Appendix B. Operator derivative relation
The differential operator L ω is linear, and we assume the boundary conditions allow for unique solutions such that the operator inverse L −1 ω exists and is well-defined.We can then write the following for the derivative with respect to τ (x): ωθ,n , explained in Appendix B. The functional derivative is thus given by δy ωθ,nm (τ, c) δτ (x) = − ∂L ω (x; τ, c) ∂τ (x) Z * ωθ,nm (x; τ, c)P ωθ,n (x; τ, c).

Figure 2 :
Figure 2: The true absorption difference profile, ∆τ = τ −τ 0 , that we wish to reconstruct.Region of interest for contrast analysis is indicated by the dashed line box where the edge of square-shaped absorption target acting as the mid-line.

Figure 5 :
Figure 5: Contrast analysis curves for 1 mm sensors.(a)-(b) show the MTF FWHM and maximum normalised contrasts for the case of 5 frequency reconstructions, while (c)-(d) show the same quantities for 1 frequency reconstructions.The solid curves correspond to PS reconstructions whereas the dashed curves correspond to the PI reconstructions.In each case there are two line pairs distinguished by the markers along the curves, corresponding to the noise-free (l) and noisy (♦) cases.

Figure 6 :
Figure 6: Contrast analysis curves for 5 mm sensors.(a)-(b) show the MTF FWHM and maximum normalised contrasts for the case of 5 frequency reconstructions, while (c)-(d) show the same quantities for 1 frequency reconstructions.The solid curves correspond to PS reconstructions whereas the dashed curves correspond to the PI reconstructions.In each case there are two line pairs distinguished by the markers along the curves, corresponding to the noise-free (l) and noisy (♦) cases.

Figure A1 :
Figure A1: Error function f (red dashed line) fitted to the ensemble average of ESFs (solid black line) for the cases of (a) d = 5 mm PS sensors with 3 angles, 1 frequency, and 1 % noise, and (b) d = 1 mm PS sensors with 24 angles, 5 frequencies, and 1 % noise.