Image reconstruction with the Heaviside equation in photoacoustic tomography accounting for dispersive acoustic media

Verena Margitta Moock; Edahí Antonio Gutiérrez-Reyes; Crescencio García-Segundo

doi:10.1117/1.JBO.23.7.076010

19 July 2018 Image reconstruction with the Heaviside equation in photoacoustic tomography accounting for dispersive acoustic media

Verena Margitta Moock, Edahí Antonio Gutiérrez-Reyes, Crescencio García-Segundo

Author Affiliations +

Journal of Biomedical Optics, Vol. 23, Issue 7, 076010 (July 2018). https://doi.org/10.1117/1.JBO.23.7.076010

Abstract

A challenging issue in photoacoustic biomedical imaging is to take into account the presence of dispersive acoustic media, since these are prone to induce amplitude attenuation and scattering of the photoacoustic frequency components. These perturbations are largely the cause for which the photoacoustic tomographic image reconstruction from projections lacks a plane-wave transport formalism. Attending this problem, we further develop an analytic formalism of the transport and its numerical implementation accounting for dispersive acoustic media. We differentiate three variations of an acoustically perturbing media. Our object of interest is a numerical description of the light absorption map of a coronal human breast image. Then, we analyze conditions for which the propagation of photoacoustic perturbations can obey the generalized Heaviside telegraph equation. In addition, we provide a study of the causality consistency of the wave propagation models. We observe transport implications due to the presence of dispersive acoustic media and derive model adjustments that include attenuation and diffusion approximations within the two-dimensional forward problem. Next, we restore the inverse problem description with the deduced perturbation components. Finally, we solve the nonlinear inverse problem with a numerical strategy for a filtered backprojection reconstruction. At a stage prior to the image reconstruction, we compensate for the effect of acoustic attenuation and diffusion to calculate the inversions of the wave perturbations located within the projections. In this way, we manage to significantly reduce reconstruction artifacts. In consequence, we prevent the use of some additional image processing of noise reduction. We demonstrate a feasible strategy on how to solve the stated nonlinear inverse problem of photoacoustic tomography accounting for dispersive acoustic media. In particular, we emphasize efforts to achieve an analytical description, and thus an algorithm is placed, for imaged sound perturbations to be cleaned from acoustic scattering in a simplified manner.

1. Introduction

The use of photoacoustic (PA) tomography in biomedical imaging has opened up interest to improve early breast cancer detection, for being a nonionizing tomography modality. The possibility of expanding the PA waveform model for acoustically perturbed media is a task of intense research. Indeed, the common plane wave model described by Xu and Wang¹ has proven to be successful in homogeneous material or at shallow depths, where dispersive acoustic media can be neglected. However, beyond the perturbation condition, the propagation model is prone to errors and is likely to introduce artifacts when reconstructing the PA source. In recent years, efforts have been focused on either reducing those artifacts that are correlated to reflections and scattering of acoustic waves or in minimizing blurring from the experimental stage.²^–⁴ Others, as Dean-Ben et al.⁵ suggested, the use of a reconstruction algorithm with a statistical detector correction of the perturbation. The authors of Refs. 6 and 7 reported how to reduce the blurring and therefore increase the contrast by taking the spatially varying sound velocity into account. Rivière et al.⁸ and Ammari⁹ in turn introduced an extended wave model whose purpose is to compensate acoustic attenuation on the measured projection information, before reconstructing the absorption distribution.

We note that so far the amendments to the underlying transport models, as implemented, still fail to account satisfactorily for recorded artifacts caused by dispersive acoustic media when working with PA sources placed within tissue at considerable depth in respect to the surface of reference for the measurements.¹⁰^–¹² These image artifacts cannot be eliminated completely, even with contrast agents;¹³ in fact, some are inherent of the actual physical properties of the medium. Along with the scenario in Ref. 14, we introduced a sound dispersion approximation as part of the PA propagation model together with attenuation considerations. This approach models the Heaviside telegraph equation in the Fourier domain. In order to deal with its inverse problems of image reconstruction related to the transport, we carried out a projection-processing strategy similar to that followed by Ammari.⁹ By following such an approach, we manage to assure that established reconstruction algorithms, such as filtered backprojection, can be applied to PA inverse problems. In the current contribution, we differentiate the acoustic perturbation among attenuation, dispersion, and a combination of both by separate operators. By these differentiations of the propagation, we are capable of simulating internal reflections and image blurring, commonly present in images produced by PA data. Furthermore, we introduce a methodology for image reconstruction in PA tomography, where we analyze, interpret, and propose a practical pathway for including fine-tuned acoustic attenuation and/or dispersion in the PA waveform model.

This paper is organized as follows: in Sec. 2, we present our proposed extension to the PA propagation model with Heaviside’s telegraph equation. This is made with the explicit purpose of accounting for dispersive acoustic media and the PA attenuation with acoustic dispersion due to the viscoelasticity properties of the transport medium; such that one can reconstruct the PA distribution source, resulting from the tissue’s absorption of energy, after being excited by laser light. We outline in Sec. 3 a strategy theory for solving the inverse problem when diffusion assumptions are made and investigate the causality condition of the models. Computer simulations of the PA inverse problem are given in Sec. 4. We present numerical solutions to distinct perturbed forward problems and demonstrate how errors, due to inappropriate projection approximation, can affect the quality of reconstruction. As part of this proposal, we also introduce a method to process the projections for perturbation adjustment. Along with our presentation, we compare backprojection results and evaluate the use of the proposed projection correction strategy. Finally, we give our conclusions and discuss briefly the implications of our model extension in Sec. 5.

2. Inverse Problems in Waveform Tomography

In mathematical terms, the PA transport is described as a system of partial differential equations, modeling the propagation of acoustic pressure distribution by the function $p : R^{n} \times R \to R_{+}$ , defined over space and time. When the initial pressure distribution $f (x) = p (x, 0)$ propagates on weakly scattering conditions (facing small objects, relative to the wavelength scale), the transport processes, in free space and with stress confinement, allow the approximation of the sound propagation as a linear integral equation $\forall x \in Ω \subset R^{n}$ and $\forall t \in [t_{1}, t_{2}] \subset R_{+}$ , which is the solution to the problem

Eq. (1)

L p (x, t) = 0,

Eq. (2)

p (x, 0) = f (x),

Eq. (3)

\partial_{t} p (x, 0) = 0 .

Given the initial and boundary conditions, the specific linear operator $L : Ω \times [t_{1}, t_{2}] \to R_{+}$ describes a linear differential equation of second order. The source distribution $f$ stands in strong physical relation to the strength of light absorption and is thus proportional with the magnitude of the PA effect.

The system of Eqs. (1) to (3) is a general description that implies distinct waveform transport analysis: this set of equations is valid and defined for a broad domain comprising acoustic, electromagnetic, optical, and seismic transport modalities. In a nonperturbing homogeneous media, the wave transport is modeled by the linear four-dimensional d’Alembert operator $□$ , which is a generalization of the Laplace operator $\nabla^{2}$ . This is so by including at once the second partial derivative in time and the constant sound wave speed $c$

Eq. (4)

L_{0} p_{0} := □ p_{0} = (\partial_{t}^{2} - c^{2} \nabla^{2}) p_{0} .

In particular for this problem, the integral operator solution of the acoustic field in dimensions 2 is given by Kirchhoff’s equation

Eq. (5)

p (x, t) = \frac{1}{2 π} \partial_{t} \int_{t_{1}}^{t_{2}} \frac{R_{Ω} [p_{0}] (x, c_{0} r)}{\sqrt{t^{2} - r^{2}}} d r,

Eq. (6)

R_{Ω} [p_{0}] (x, c_{0} r) = 4 r \int_{0}^{r} \frac{p (x, t / c_{0})}{\sqrt{r^{2} - t^{2}}} d t .

However, for acoustic dispersive media, the above plane wave propagation is erroneous, since its lack of consideration of the diffusion and attenuation processes, both physically intrinsic conditions of the general phenomena.¹^,¹⁵ When attenuation is present, the waveform operator expresses the Helmholtz equation as applied in ultrasound imaging, and even PA imaging with a complex wave number.¹⁶ The corresponding mathematical expression is given as

Eq. (7)

L_{a} p_{a} := □ p_{a} + a * p_{a} .

When sound dispersion occurs as in optical tomography,¹⁷ the Boltzmann equation applies

Eq. (8)

L_{d} p_{d} := □ p_{d} + d * \partial_{t} p_{d} .

In Eqs. (7) and (8), $a$ and $d$ are appropriate weight functions over $C$ and $*$ denotes the convolution between two functions. In more complex cases, attenuation and dispersion are simultaneously present. This situation is common in PA imaging of deeper biological tissues.¹⁸ The condition leads to the Heaviside telegraph equation

Eq. (9)

L_{a, d} p_{a, d} := □ p_{a, d} + d * \partial_{t} p_{a, d} + a * p_{a, d} .

This model has been tested by Arridge¹⁹ in the field of optical tomography. Moock et al.¹⁴ explored this scenario, with the aim of modeling the PA transport for dispersive acoustic media, which is a more realistic representation of the breast tissue. The study prompted to disclose different wave conditions arising from the inclusion of acoustic perturbations.

In waveform tomography, it is assumed that for any source $f$ there exists a projection $g$ dependent on the detector position that corresponds to what acoustic sensors register at the boundary of the observed closed and bounded region $Ω$ ; over the fixed time interval $[t_{1}, t_{2}]$ . The inverse problem of this waveform tomography, in two steps, states:

a. Given the projections at the boundary $p (y, t) = g (t)$ , with $y \in \partial Ω$ , and $t \in [t_{1}, t_{2}]$ .
b. Find the initial sound pressure distribution $f (x)$ , whose projection is then, $g (t)$ .

When it comes to solve the inverse problem, it is important to observe that although the transport obeys a propagation model of linear equations, the related inverse problem may be nonlinear. In fact, in the case of PA, the inversion of the operator $L$ involves the pressure distribution $p$ and other additional weight functions $a, d$ , which cannot be reconstructed at once, e.g., by backprojecting the measured data $g (t)$ .

3. Solution of the Dispersive Photoacoustic Inverse Problems

The inversion of the PA perturbations is the central task of the considered dispersive inverse problems modeled in Sec. 2. For biological tissues such an endeavor becomes difficult since linear operators $L_{a}$ , $L_{d}$ , and $L_{a, d}$ impose a nonlinear inverse problem. Ref. 14 followed the practical strategy of numerically solving the PA inverse problem modeled by Eq. (9) for an acoustically dispersive medium. This method encourages signal preprocessing including some boundary measurements. As a result, the impact of sound perturbations within the media tends to disappear. In the current contribution, we compare the three propagation models Eqs. (7)–(9) for image reconstruction on synthetical two-dimensional (2-D) data sets.

When applying the strategy from Ref. 14 to all waveform conditions, the formalism requires the description of Eqs. (4)–(9) in the Fourier domain, and consequently with the pressure distribution expressed by a frequency-dependent function, $\hat{p}$ . Then, the corresponding operators are

Eq. (10)

(c^{2} \nabla^{2} + ω^{2}) {\hat{p}}_{0} (x, ω) = \frac{\hat{ı} ω}{\sqrt{2 π}} f (x),

Eq. (11)

[c^{2} \nabla^{2} + ω^{2} - \hat{a} (ω)] {\hat{p}}_{a} (x, ω) = \frac{\hat{ı} ω}{\sqrt{2 π}} f (x),

Eq. (12)

[c^{2} \nabla^{2} + ω^{2} - \hat{ı} ω \hat{d} (ω)] {\hat{p}}_{d} (x, ω) = \frac{\hat{ı} ω}{\sqrt{2 π}} f (x),

Eq. (13)

[c^{2} \nabla^{2} + ω^{2} - \hat{ı} ω \hat{d} (ω) - \hat{a} (ω)] {\hat{p}}_{a, d} (x, ω) = \frac{\hat{ı} ω}{\sqrt{2 π}} f (x) .

Since in a homogeneous media, there is no random field, the pressure field is the coherent part to the pressure distribution ${\hat{p}}_{0}$ in Eq. (10) and ${\hat{p}}_{a}$ in Eq. (11). Instead, although Eq. (12) is similar to the electromagnetic diffusion in a conducting medium, the field ${\hat{p}}_{d}$ should be interpreted as the energy density of the pressure field and aims to account for the diffused part of the acoustic field. In Eq. (13), the field ${\hat{p}}_{a, d}$ should be interpreted as the energy density of the pressure field, covering the attenuated and diffused parts of the acoustic field. The coherent part is ignored because we suppose the scatters concentration to be dense enough as to inhibit it.²⁰

For each case, the source term $f (x)$ can be interpreted as the initial pressure, but if we allow a frequency dependence $f (x, ω)$ , then it can be proportional to the time derivative of the temperature field if heat transport is considered, or proportional to the electromagnetic density energy if no heat transport is considered.⁴^,²¹^,²²

In this formalism, one can derive an expression between the homogeneous and the boundary measurements of dispersive acoustic media, and then obtain their respective wave numbers, $k_{0}, k_{a}, k_{d}$ , and $k_{a, d}$ . The homogeneous case would produce $k_{0} = \frac{ω}{c_{0}}$ , the wave number of a perfect monochromatic wave, also known as the dispersion relation. In it, $c_{0}$ is the constant sound speed of the given homogeneous medium. However, when the field is nonmonochromatic and propagates with conditions of acoustic attenuation with dispersion and scattering, the complete dispersion relations for the propagation of the field as described from Eqs. (10)–(13) is

Eq. (14)

k = \frac{| ω |}{c_{o}} \sqrt{1 - \hat{ı} \frac{\hat{d} (ω)}{ω} - \frac{\hat{a} (ω)}{ω^{2}}},

all different cases are properly covered using the proper constant values, as noted in Ref. 14. Here we need to make a compromise between having a logic consistency and stable numerical results. For the calculation of the propagating field, we solve the problem classically for Eqs. (10)–(13), by means of using their corresponding integral solutions to state the inverse problem; instead, for compensating the attenuation we introduce as novelty, the use of empirical attenuation relations.

3.1.

Sound Velocity as Function of the Frequency

Among the several physical properties that characterize a colloid-like condensed medium, such as the biological tissue, one can account for its density and its viscoelastic properties. In terms of PA interpretation, if the medium is homogeneous, it can be characterized by its attenuation and acoustic dispersion. However, when the medium is acoustically dispersive, and simultaneously exhibit scattering, both properties prevent the sound velocity to remain constant over the transport process. Therefore, in either case we can assume a physical sound velocity distribution that is characteristic of the given medium. Although one can get to know the characteristic PA acoustic spectral distribution, a priori we cannot predict the spatial evolution of such spectral distribution over the transport process. Therefore, for simplicity we assume that the simplest form of the dispersion relation is given as

Eq. (15)

k (ω) = \frac{ω}{c_{p} (ω)},

where

c_{p} (ω)

is the phase velocity that depends on the frequency. As noted elsewhere in the literature in Refs. 5, 8, 9, 14, 20, 23, and 24, the initial PA perturbation in the frequency domain is given by an initial wave packet centered around a certain frequency. However, when the perturbation propagates under conditions of either, attenuation with acoustic dispersion, or changes of acoustic impedance, or both, then the frequency distribution gets attenuated and the wave packet gets dispersed. Thus, depending on the length of propagation and on the density and magnitude of impedance changes, the velocity of propagation of the wave packet would be dominated by a group velocity, which in turn is given by the derivative of the wave vector with respect to the change in frequency. Therefore, in the classification noted above, this dependency is implicit through the respective wave number. Even so, in general, one cannot get to know a prior which one of the noted perturbations would cause the frequency changes. However, for the theoretical representation purpose, these can be treated as dispersive acoustic media. Thus, one can model conditions that inhibit the PA perturbation of keeping a constant velocity and to be monochromatic, and yet fulfill the wave equation and introduce means for improving the image reconstruction process, as described later in Secs. 3.2 and 4.

As reported in Refs. 14, 16, 23, and 24, the change of depth, density (say, viscoelasticity), or acoustic impedance would force the signal propagation to change its velocity distribution (group velocity). In Ref. 14, we verified experimentally the changes in the frequency spectrum, which in the temporal domain are expressed as temporal stretching or broadening, as a function of the boundary conditions and physical transport properties.

Our stated proposition derived from these observations is that the wave vector is frequency dependent, and in consequence, the sound velocity cannot be constant. Since we cannot get to know the actual propagation trajectory, including specific impedance changes, we suggest that through a global modulation, via the frequency-dependent sound velocity, one can summarize the actual wave vector and the effective sound velocity. Therefore, we suggest to consider a group-velocity-like, that is, modulating the wave packet that gives shape to the PA mechanical perturbation. An analogous case to this phenomena is related to the propagation of an optical field through a random medium, such as a colloid, where the wave number turns out to be a complex function,²⁵ which is analogous to this displayed here in the Eq. (10). In Refs. 14 and 23, we displayed and discussed how is that the sound speed keeps drifting the bandwidth and the peak value, rather than keeping a constant speed distribution. Hitherto the inverse problem and thus the conventional PA tomography reconstruction does not incorporate this change in sound velocities as part of the solution strategies. Therefore, one of the central motivations of the current contribution is to place a formal discussion of the subject, and then work out a strategy to dealt with this issue in a practical manner, specifically related to the image reconstruction in PA tomography. To the best of our knowledge, we assess that several authors have introduced methodologies using additional experimental features to include reference signals or else treating these in a postreconstruction stage (for the PA signals or the PA tomography images), to compensate specific image reconstruction defects. Yet, none of these methods introduce or treat the problem in the way we dealt with it. The interested reader can review for instance.²^–⁸^,¹⁶

In this work, we propose a practical projection approximation to treat and compensate dispersive acoustic media causing either attenuation, or sound dispersion, or both, directly within the recorded PA data. The homogeneous transport gives the recorded data $g_{0}$ as the projection of ${\hat{f}}_{0}$ , with wave number $k_{0}$ . As described by Gutierrez-Reyes et al.,²⁴ by means of using the appropriate sensing process, the signal read out would not alter the projection characteristics in a meaningful way. In the presence of dispersive acoustic media, the PA perturbations would be distorted due to changes of the acoustic impedance (acoustic scattering) and attenuation with dispersion due to transport length in a viscoelastic medium. Both contributions impose changes on the frequency distribution and, in consequence, on the characteristic velocity distribution, depending on the acoustic dispersive properties of the media. This behavior is apparent from the frequency-dependent wave number $k_{a} (ω), k_{d} (ω)$ , or $k_{a, d} (ω)$ .¹⁴ Thereafter, when referred generically to the distinct type of perturbations, we will use the subindex $q = {a, d, a d}$ , respectively. Otherwise, each distinctive index would be used for naming each perturbation type. Consequently, we treat the signal alteration as a projection transformation, caused by the transport perturbation $g_{q}$ , accordingly with the $q -$ perturbation type. We define $g_{q} = T_{q} g_{0}$ , where

Eq. (16)

T_{q} g_{0} (t) = \frac{1}{2 π} \int_{R} \frac{ω}{c_{0} k_{q} (ω)} e^{\hat{ı} ω t} \int_{t_{1}}^{t_{2}} g_{0} (t^{'}) e^{- \hat{ı} c_{0} k_{q} (ω) t^{'}} d t^{'} d ω

represents how to transform the undistorted projection

g_{0}

, to become the perturbed projection

g_{q}

. As noted before, the

q

-subindex runs as

q = {a, d, a d}

, depending on the perturbation that is considered: attenuation, sound dispersion, or both, respectively. Since

c_{0}

is constant, without loss of generality we can assume

c_{0} = 1

. We notice that Eq. (16) is a consequence that in the Fourier space

g_{0}

and

g_{q}

are related as

Eq. (17)

{\hat{g}}_{0} [x, k_{q} (ω)] = \frac{k_{q} (ω)}{ω} {\hat{g}}_{q} (x, ω) .

The first integral implicitly introduces the attenuation and diffusion to the clean measure $g_{0}$ by evaluating its Fourier transform at $k_{q} (ω)$ , while the second integral (left) is the inverse Fourier transform of Eq. (17). In practice, we know the pressure $g_{q}$ at the boundary and want to recover $g_{0}$ also at the boundary. This is done by inverting the transformation Eq. (16). Further, the three variants of the dispersive acoustic media, represented by Eq. (16), require some knowledge of the function $k (ω) = k_{q} (ω)$ . This function was discussed in Ref. 14, where it was disclosed that, depending on the perturbation type, it can be a real or a complex number. Respectively, the wave number in the three study cases of acoustical perturbation we discuss becomes

Eq. (18)

k (ω) \leftarrow {\begin{cases} k_{a} (ω) = \frac{ω}{c_{0}} + \frac{1}{2} \hat{ı} \bar{a} ω^{2}, \\ k_{d} (ω) = \frac{ω}{c_{0}} - \frac{1}{2} \bar{d} ω^{3}, \\ k_{a, d} (ω) = \frac{ω}{c_{0}} + \frac{1}{2} (\hat{ı} \bar{a} ω^{2} - \bar{d} ω^{3}) . \end{cases}

Here, $\overline{a}$ and $\overline{d}$ represent appropriate weighting factors, mimicking perturbations in the tested transport media. The experimental model conditions in Ref. 14 suggest equal constant values for $\overline{a}$ and $\overline{d}$ . Here, we simulate a PA scenario at real breast scale (the diameter of the average female human breast measures 25 cm). Accordingly, we set the perturbation weights $\overline{a} = - 0.1$ and $\overline{d} = 0.1$ . Moreover, these parameters are close to the model settings of Ref. 9 and keep all calculations stable.

3.2.

Causality Condition

Here, we investigate the causality of the former models. We can assume that we have a dispersive medium where the wave vector is $k (ω) = \frac{ω}{c_{p} (ω)}$ , as previously noted in Eq. (15), where we allow the velocity to depend on the frequency, and this being a causal response implies the existence of an imaginary velocity related to the real part by the Krammers–Kronig relations that can be written as a Hilbert transform. One needs to take care of whether the Krammers–Kronig relations apply or not to the wave vector. In general, the wave vector depends on the incidence angle and boundary conditions, which are not intrinsic material properties. Therefore, rather we consider it to be a self consistent value for a certain wave equation with dispersion. To verify this assessment we discuss three cases referring to the wave numbers of Eq. (18). The first one is already well known in the literature⁹^,²⁶ and the second has been partially discussed in Refs. 17 and 19, while the third one relates to the Heaviside equation and that is central for the current contribution. These cases are as follows.

3.2.1.

Attenuation with acoustic dispersion

Here the wave vector magnitude is defined as

Eq. (19)

k (ω) = \frac{ω}{c_{p} (ω)} = \frac{ω}{c_{0}} + \frac{\bar{a}}{2} \hat{ı} w^{2} .

We know from the properties of the Hilbert transform that for a function of complex variable such that the imaginary part is the Hilbert transform of the real part, it should be an analytic function (as noted in Sec. 15.1.2 in Ref. 27). Minding Eq. (19), complex analytic function for the phase of the sound speed is given as

Eq. (20)

c_{p} (z) = \frac{1}{\frac{1}{c_{0}} + \frac{\bar{a}}{2} \hat{ı} z}

and is an analytic function with

z = ω + \hat{ı} τ

. When

τ = 0

we find that

Eq. (21)

c_{p} (ω) = \frac{1}{\frac{1}{c_{0}^{2}} + \frac{{\bar{a}}^{2}}{4} ω^{2}} (\frac{1}{c_{0}} - \frac{\bar{a}}{2} \hat{ı} ω) .

Therefore, the quadratic model describes a material whose real part depends on the frequency, which is represented by a Cauchy distribution; turning out then that the parameter $\overline{a}$ is related to the width of the distribution, in which the imaginary part is its Hilbert transform.

3.2.2.

Acoustic diffusivity

This time, let us assume that the velocity is real, although diffusive as a result of scattering due to acoustic impedance changes. It places conditions for the PA-induced mechanical perturbation that prevents it from being monochromatic and in consequence, an error in causality is introduced as

Eq. (22)

k (ω) = \frac{ω}{c_{p} (ω)} = \frac{ω}{c_{0}} - \frac{\bar{d}}{2} w^{3} .

In these conditions we can choose its real part

Eq. (23)

Re [c_{p} (ω)] = \frac{1}{\frac{1}{c_{0}} - \frac{1}{2} \bar{d} ω^{2}},

whose Hilbert transform can be found using a partial fraction expansion. Furthermore, after using the result Eq. (32) in Table 15.2 from Ref. 27, it turns out that if

\overline{d} > 0

, then the imaginary part is

Eq. (24)

Im [c_{p} (ω)] = \frac{π \sqrt{c_{0}}}{\sqrt{2 \bar{d}}} [δ (ω - \sqrt{\frac{2}{c_{0} \bar{d}}}) - δ (ω + \sqrt{\frac{2}{c_{0} \bar{d}}})] .

However, if instead $\overline{d} < 0$ , then its Hilbert transform gives

Eq. (25)

Im [c_{p} (ω)] = - \frac{\sqrt{\frac{c_{0} d}{2}} w}{\frac{1}{c_{0}} + \frac{1}{2} d ω^{2}} .

Therefore, without an imaginary part, this model fails to fulfill Krammers–Kronig and hence the causality condition. This case resembles an analogy to the well-known anomalous optical dispersion. So far, up to the best of our knowledge, we recognize it as yet an open problem in PA interpretation.

3.2.3.

Attenuation with acoustic dispersion and diffusive scattering

This case is central for the present discussion. Its wave vector’s magnitude is set to

Eq. (26)

k (ω) = \frac{ω}{c_{p} (ω)} = \frac{ω}{c_{0}} + \frac{1}{2} (\hat{ı} \bar{a} ω^{2} - \bar{d} w^{3}) .

Here we can choose

Eq. (27)

c_{p} (z) = \frac{1}{\frac{1}{c_{0}} + \frac{1}{2} (\hat{ı} \bar{a} z - \bar{d} z^{2})}

with

z = ω + \hat{ı} τ

, which separating the real part and imaginary part and taking

τ = 0

we get that

Eq. (28)

c_{p} (ω) = \frac{\frac{1}{c_{0}} - \frac{1}{2} \bar{d} w^{2} - \hat{ı} \frac{1}{2} \bar{a} w}{{(\frac{1}{c_{0}} - \frac{1}{2} \bar{d} w^{2})}^{2} + \frac{1}{4} {\bar{a}}^{2} ω^{2}} .

Thus, causality prevails. Despite this fact, it becomes apparent that the imaginary part in the velocity spectra model is required to describe causality with physical meaning. Instead, if the velocity is a real domain parameter (depending on the frequency), then an error is introduced that fails the causality condition. Noticeably, one can get an explicit expression for the frequency-dependent velocity distribution. In further sophistication, one can expect to properly describe the coefficients $\overline{a}$ and $\overline{d}$ , and find a proper analytic description in terms of physical properties. However, such a task is of a nature beyond the current scope and would be treated in a separate analysis.

By the analysis of the three cases, the constant values $\overline{a} = - 0.1$ and $\overline{d} = 0.1$ do not affect the general argument of causality. These values are only for the purpose of proof of concept and testing the performance of the interpretation and reconstruction model. We agree that a much more detailed description and precise range of values for $\overline{a}$ and $\overline{d}$ is desired. However, the task by its nature demand deep mathematical analysis to obtain the conditions, and range of values and solutions that are stable and have a correct physical meaning, and at once also satisfy the requirements for achieving satisfactory image reconstructions.

4. Photoacoustic Simulations

4.1.

Forward Problem

In order to demonstrate how the modeled sound perturbations affect the PA reconstruction, we carried out a simulation of the 2-D propagation through the phantom in Fig. 1. This phantom resembles the anthropomorphic features of a cross section of a female human breast according to Ref. 28. It is discretized to offer a simple visualization of different constant geometries (i.e., indicated by manual annotations of medical experts) and is here described as a collection of several overlapped circles and one star displayed by gray scales at distinct values from $[0,1] \cap Q_{+}$ . The data set was generated with a stable and reliable open source vector graphics editor in such a way that it captures the geometrical properties of the characteristic biological breast components. Based on Ref. 28, we assign normalized optical absorption coefficients to the different regions values in the range [0.5, 1] to specify sections of major absorption (dark gray picture elements), as it is the case of cancerous and glandular tissue. On the other hand, values in the range [0, 0.5] are to specify sections of lesser absorption (brighter picture elements), such as fat and fibroadenoma. According to the PA effect, the light absorption map resembles the initial pressure distribution. The soft tissues that we include in the phantom appear typically in a common mammograph and similar values are expected for PA tomography due to familiar energy absorption characteristics; as an example, see the results shown in Ref. 29. For testing our algorithms, we used a digitized version of the phantom using a $400 \times 400$ image, which allows enough resolution for the details of the specimen.

Fig. 1

A $400 \times 400$ digitization $f$ of the phantom $f$ representing structures in the coronal view of a breast by geometrical shapes. The light gray circles (a), (b), and (c) represent bounded and unbounded fibroadenomas and the dark circles (d) and (e) represent small cancerous tumors; a larger tumor is here depicted by a star (f). The expected absorption coefficients in the anthropomorphic features are quantified with a value arbitrarily assigned from the interval $[0,1] \cap Q_{+}$ , while the expected absorption coefficients in the breast glandular tissue (g) are assigned to a value equal to 0.8 and in fat (h) to 0.2 (in arbitrary units).

Once we decided on the corresponding sound attenuation and dispersion weights $\overline{a}$ and $\overline{d}$ for a specific application, and on the complex wave number $k (ω)$ in Eq. (18), the perturbed forward solutions $g_{a}, g_{d}$ , and $g_{a, d}$ can be easily derived by the homogeneous solution $g_{0}$ according to the mathematical expressions in A1-A3 of Sec. 3. For our experiments, we use 180 parallel projections acquired every degree with 400 lines each. The resulting projection integrations are approximated by $g_{0} \approx R f$ . The perturbed projection data are derived by the transformation $T_{q} g_{0}$ as outlined in Eq. (16); i.e., $g_{a} \approx T_{a} g_{0}$ in the case of acoustic attenuation, $g_{d} \approx T_{d} g_{0}$ in the case of sound dispersion, and $g_{a, d} \approx T_{a, d} g_{0}$ in the case of both perturbation factors.

4.2.

Classical Filtered Backprojection Results

The digital reconstruction ${f^{'}}_{q}$ produced by the classical filtered backprojection algorithm is obtained by applying inverse Radon transform with a Ram–Lak filter to the calculated projections $g_{q}$ ; i.e., ${f^{'}}_{q} = R^{- 1} g_{q}$ . The projections are fast-Fourier transformed and multiplied with the ramp filter in the frequency domain. After backprojecting every set of the differently perturbed data, one can simulate the known artifacts: internal reflections⁵ and image blurring,⁶^,⁷ see Fig. 2. From the visual inspection, we can give the following observations: while the result in Fig. 2(a) is an almost error-free image of the original pressure distribution map $f$ , it poorly represents a realistic reconstruction of a breast by PA imaging in comparison to realistic reconstructions of Ref. 30. Backprojection of $L_{a}$ results in Fig. 2(b) bringing out a blurring effect increasing toward the upper half of the image, generated in accordance to our attenuation approximation. On the other hand, backprojection of $L_{d}$ results in Fig. 2(c) producing internal reflections, directed toward one side, as determined by the transport propagation direction. As expected from the telegraph equation, Fig. 2(d) describes a reconstruction with both phenomena, blurring and less visible reflections. Besides one can assess that the blurring component turns out to be dominant over the reflections. A task oriented figure-of-merit (FOM) is suitable for a qualitative evaluation of algorithm performance. FOMs provide different measures of picture distances between a reconstruction and a phantom. Next, we recall two of the most known standard criteria for qualitative assessment of the reconstruction.³¹^,³²

Fig. 2

Digital reconstructions ${f^{'}}_{q}$ from the projection data $g_{q}$ : (a) no transport-dependent artifacts are visible in ${f^{'}}_{0}$ when acoustic perturbations are absent (in the case of $g_{0}$ ); (b) ${f^{'}}_{a}$ shows image blurring; (c) ${f^{'}}_{d}$ shows internal reflections; and (d) ${f^{'}}_{a, d}$ shows the combination of reflections and blurring in the same reconstructed image.

4.2.1.

Normalized root mean squared distance

This criterion emphasizes the importance of large errors throughout all the $J$ pixels of the discretized images

Eq. (29)

d_{RMS} (f, f^{'}) := \sqrt{\frac{\sum_{j = 1}^{J} {(f_{j} - {f^{'}}_{j})}^{2}}{\sum_{j = 1}^{J} {(f_{j} - \bar{f})}^{2}}},

where

f_{j}

and

{f^{'}}_{j}

indicate the

j

’th element of

f

and

f^{'}

, respectively. A small value of

d_{RMS} (f, f^{'})

indicates a small reconstruction error.

4.2.2.

Structural accuracy (STR)

Let ${\bar{f}}^{(m)}$ and ${\bar{f}}^{' (m)}$ be the average pixel value of $f$ and $f^{'}$ , respectively, of those pixels whose centers are within structure $m$ from the total of $M$ structures. Then, we define the structural accuracy as

Eq. (30)

d_{STR} (f, f^{'}) := - \frac{1}{M} \sum_{m = 1}^{M} | {\bar{f}}^{' (m)} - {\bar{f}}^{(m)} | .

For the current case, we compare a total number of six structures ( $M = 6$ ) as can be distinguished in Fig. 3.

Fig. 3

Six objects of the phantom of Fig. 1 for a structural accuracy evaluation.

In Table 1, we show the numerical comparison of the perturbed backprojection displayed in Fig. 2. The numerical values are obtained from the above-mentioned FOMs. We notice that the acoustic dispersion deteriorates mainly the structural accuracy and increases the normalized root mean squared (RMS) distance between the digitized phantom and its reconstruction. Interestingly, in this specific reconstruction example, the modeled acoustic perturbations are not linearly dependent. We consider that this case represents a realistic situation because attenuation and dispersion have physically distinct origins and implications.

Table 1

Mean distance and structural performance of the perturbed noise-free backprojection.

Transport	dRMS (f,f′q)	dSTR (f,f′q)
$L_{0}$	0.1770	$- 0.0160$
$L_{a}$	0.2990	$- 0.0623$
$L_{d}$	0.4153	$- 0.0546$
$L_{a, d}$	0.2944	$- 0.0588$

4.3.

Corrected Photoacoustic Projection Data and Improvement of the Image Reconstruction

The solution strategies elaborated in Sec. 3 give a hint about how to obtain an estimation of the transport inversion. For further simulations in a computer, we need to discretize the transport $T$ in propositions A1-A3 with the complex wave-number description in Eq. (18); however, its analytical inversion cannot be guaranteed since the linear system is ill-conditioned because of the perturbed description. An acceptable discretization $T$ for the three propositions of A1-A3 can be derived by numerical analysis. One such discretization is the singular value decomposition (SVD), which approximates $T$ by the matrix product ${UWV}^{t}$ where $U$ and $V$ are real or complex unitary matrices (whose transpose is represented by $V^{t}$ ), and $W$ is a rectangular diagonal matrix. For a positive real number $ϵ > 0$ , we can obtain an approximation of the inverse of $T$ according to

Eq. (31)

T_{1, ϵ}^{- 1} g = \sum_{j = 1}^{J} \frac{w_{j}}{w_{j}^{2} + ϵ^{2}} ⟨ g, V_{j} ⟩ U_{j}^{t},

with

U_{j}

and

V_{j}^{t}

representing the

j

’th row of matrices

U

and

V^{t}

, respectively, and all

w_{j}

, for all pixel indices

1 \leq j \leq J

, are nonnegative real singular values (i.e., the square roots of the

J

eigenvalues). As a follow-up study of the simulations, the inversion of the discretized transport description of propositions A1-A3 to better approximate projection data was applied in Ref. 14. In fact, this strategy has also been used by other authors, such as in Ref. 9. In respect to the SVD algorithm, when used for this purpose, it is known to become unstable when

ϵ

tends to zero. Hence, in order to optimize the inversion, we applied a simple line search strategy, which converged at

ε = 20

.

We want to demonstrate that the corrected projection data are in agreement with the PA transport model, i.e., $T_{a}^{- 1} g_{a}$ instead of $g_{a}$ when attenuation is present, $T_{d}^{- 1} g_{d}$ instead of $g_{d}$ when dispersion is present, and $T_{a, d}^{- 1} g_{a, d}$ instead of $g_{a, d}$ in the case of both modeled acoustic perturbations and it lead us to improved backprojection results. Since a qualitative evaluation is restricted to a computer simulation, as performed here, we tried to incorporate a realistic error description into our mathematical model.

The signal-to-noise analysis in Ref. 33 provides a noise description of biomedical PA measurements and indicates how to integrate this noise into our simulation. We have tested different piezoelectric sensors and measured their ambient noise. Consequently, we added to our projections additive Gaussian white noise with mean $μ = 50 \times 10^{- 4}$ and standard deviation $σ_{1} = 2 \times 10^{- 4}$ . The resulting signal-to-noise ratio (SNR), as defined in Ref. 34, is equivalent to $μ / σ_{1} = 25 \approx 14 dB$ , which is similar to the case study of Ref. 33. In addition, a realistic scenario is to assume that light scattering in deeper biological tissues produces a negative impact on the PA SNR.²² Therefore, we looked at the scenario when the SNR is reduced to 25%; thus, we selected $σ_{2} = 8 \times 10^{- 4}$ , which results in $μ / σ_{2} = 6.25 = 10 \log_{10} 6.25 dB \approx 8 dB$ .

In Fig. 4, we present various backprojection reconstructions obtained from inverting the simulated forward results, and after following the strategies A1-A3. In the left column, we show the results of backprojecting the perturbed data while in the right column we show the results after appropriate adjustments that compensate for the perturbation of the modeled dispersive acoustic media. A visual inspection of these results suggests that our methodology was able to minimize both image blurring and internal reflections; the latter problem seems to be particularly well resolved in this computer simulation.

Fig. 4

Digital image reconstruction $f^{'} = R^{- 1} T_{q}^{- 1} g_{q}$ : on the left side are the results of the common filtered backprojection in the presence of additive Gaussian white noise with an SNR of 14 dB. On the right side are the backprojection after adjustments to th projection data compensating the perturbation of dispersive acoustic media. The correspondence with Fig. 2 is placed as follow: (a) and (b) correspond to blurring from attenuation, (c) and (d) to acoustic reflections, and (e) and (f) include both attenuation and reflections.

For a more quantitative way of evaluating the image reconstructions, we compared the plots of pixel densities for a representative column in the phantom and the reconstruction, as shown in Figs. 5 Fig. 6–7. In all studied cases, the reconstruction error, defined as the sum of all differences between the reconstruction using corrected PA projection data (black dash-dot curve) and the ground truth (blue continuous curve), has been considerably reduced, in comparison with the standard delay-and-sum algorithm (red-dashed curve). When the attenuation perturbs the transport, the discrepancy between both backprojections and phantom is less noticeable. The backprojections of the sound dispersion in turn bring out a strong reduction of the error oscillations. In the presence of noise, we perform qualitative image reconstruction analysis taking into account the previously mentioned FOMs. The results of this analysis are listed in Table 2 (average values), where the backprojection and the perturbation adjusted reconstructions are contrasted. It is apparent that the application of transport inversion led to success in almost every transport study; the exception is the case of $L_{a}$ with the poor SNR of 8 dB. In particular, when the dispersive acoustic media provokes internal reflections due to acoustic perturbations, the proposed projection adjustment is highly recommended. When the SNR is high, the sound dispersion correction achieved a noticeable error minimization (especially in the center of the reconstruction) with respect to both quality measures $d_{RMS}$ and $d_{STR}$ . In turn, the modeled Gaussian white noise has a mean value close to zero. This correction became less noticeable as the SNR diminished. From these results, we may infer that an inversion of the noisy projection information under the specific circumstances does not adversely affect the image reconstruction. On the contrary, the projection processing strategy provides a fundamental improvement on image reconstructions in the presence of acoustic perturbations. Here, we identified the model inconsistencies and corrected the projection data at the stage prior to reconstruction.

Fig. 5

Transport assumption $T_{a}$ with different SNR estimation: (a) 14 dB and (b) 8 dB. The continuous curve (blue) is the ground truth, the dashed curve (red) is the delayed and sum case, and the dash-dot curve (black) is the backprojection with signal processing.

Fig. 6

Transport assumption $T_{d}$ with different SNR estimation: (a) 14 dB and (b) 8 dB. The continuous curve (blue) is the ground truth, the dashed curve (red) is the delayed and sum case, and the dash-dot curve (black) is the backprojection with signal processing.

Fig. 7

Transport assumption $T_{a, d}$ with different SNR estimation: (a) 14 dB and (b) 8 dB. The continuous curve (blue) is the ground truth, the dashed curve (red) is the delayed and sum case, and the dash-dot curve (black) is the backprojection with signal processing.

Table 2

Mean distance and structural performance of the perturbed noisy backprojection.

Transport	dRMS(f,·)		dSTR(f,·)
Transport	R−1g	R−1T−1g	R−1g	R−1T−1g
$SNR = 14 dB$
$L_{a}$	0.2874	0.2100	$- 0.0544$	$- 0.0270$
$L_{d}$	0.4220	0.1906	$- 0.0561$	$- 0.0181$
$L_{a, d}$	0.2863	0.2021	$- 0.0532$	$- 0.0235$
$SNR = 8 dB$
$L_{a}$	0.2822	0.2860	$- 0.0475$	$- 0.0430$
$L_{d}$	0.4174	0.2374	$- 0.0541$	$- 0.0281$
$L_{a, d}$	0.2852	0.2807	$- 0.0518$	$- 0.0380$

5. Discussion and Conclusion

In this paper, we have extended the previous work,¹⁴ detailing concepts hinted at and however left aside. Among the most noticeable contributions we cover now is the introduction of the phase velocity [see Eq. (28)], with explicit frequency dependence. As a requirement, we also covered the causality conditions that apply to the different wave functions employed to describe the PA signal propagation. Also as extension for the theoretical interpretation, we suggest an integral transformation as a tool for correcting distortions induced over the initial PA source; which can also be understood as a density of energy, with temporary influence from a locally induced temperature field.²¹ In addition, we clarify the model description, as a diffusive model of the pressure field in contrast with coherent model description. This consideration encourages an open discussion for a better interpretation of the range of parameters involved in the propagation interpretations. The improvement is therefore more on the understanding of the meaning and application of the concepts, rather than a fundamental change in the state-of-the-art PA tomography technique, which is already highly sophisticated.

By now, we have demonstrated that the PA transport in dispersive acoustic media can be modeled with the Heaviside telegraph equation, from which can be introduced a consistent criterion to realize PA image reconstruction accounting for dispersive conditions of a given acoustic media, in a simplified manner. Among the implications we show that if the projection data are identified as the Radon transform of the perturbed source, it leads to bad approximations of the initial pressure [or absorption distribution, or time derivative of the temperature, according to the interpretation we assign to $f (x, ω)$ ]. Moreover, by a computational simulation we demonstrated how this misinterpretation affects the reconstruction finesse. Yet this apparent inconvenience brought us to introduce a method for achieving significant improvement in all study cases. This is made after using a quality and performance evaluation of the backprojection algorithm; i.e., by preprocessing the acquired data we manage to compensate the lack of information content, which is a consequence of the natural loses over the transport process. The general assessment is that the improvements we display here are introduced without demanding further computing efforts, compared with existing methods.

The findings presented here show the applicability of common algorithms for image reconstruction from projections to PA imaging. The general requirement is to ensure that the acquired data comply with an appropriate perturbed transport model, as suggested in Ref. 9. Thus, our work encourages the usage of recent improvements on reconstruction algorithms to produce much improved results. Consistent with most PA transport considerations, our model extension unifies actual theoretical and experimental findings of PA imaging theory and provides tools for interpreting noise and other sources of image artifacts. Along the way, we place a classification of the waveform transport, which are inspired by PA methods that potentially will benefit biomedical imaging.

Another important aspect of the current contribution is the introduction of the phase and group velocity concepts; we define an intrinsic sound velocity depending on the frequency, that is compatible with the causality condition (possibly a complex number). As a preliminary step to introduce these concepts, we choose the parameters according to good numerical results supported by previous reports.⁹ In reality, the parameter $\bar{a}$ should be a very small number, $\sim 6.4 \times 10^{- 14} s^{2} m^{- 1}$ for water at 20°C as in Ref. 35 and positive, and it is well understood that it arises from thermal conductivity and viscosity.²⁶^,³⁶ The other parameter $\overline{d}$ is related to the diffusion coefficient; and therefore, somehow is linked to the mean free path between scatterers;²⁰ so far we must consider it just as an empirical model, because we are not using the exact dispersion relation and therefore the meaning has not yet been clarified. Nevertheless it is noticeable that, when it is present, the meaning of $\overline{a}$ changes. As we see, when $\overline{a} = 0$ and $\overline{d} \neq 0$ , the model description fails to be causal. Indeed this subject requires more research, aiming among other things, to find out these conditions for which the current predictions match the physical properties in a close manner. By now we posed it as an intermediate step to model a diffusive media, keeping in mind that this model stage can be verified experimentally, and their proper range of values determined. We also want to emphasize that derived from the present results, the approach for compensating diffusion and attenuation using the complete dispersion relation is a topic that deserves attention. Derived from such considerations we managed to achieve further stable improvements on the task of image reconstruction; in terms of quantity and quality of information, as a complement to the existing theory. As noted within the contribution, there are several other aspects to be furthered and thus assess, which other outcomes can be gained.

Indeed it is clear that our proposed methodology can be extended to other modalities of waveform tomography by correspondingly estimating the model parameters for the specific inverse problem. Once the appropriate values are selected, it is important to adjust the acquired PA data to the perturbed transport model. We are aware that our proposal for PA tomography represents a partial step, and further work to characterize the noise present in the data is yet required in terms of different kinds of sensing methods.²⁴^,³⁷ Advances in this problem with different perspective in mind may also be found on the field of inverse media or obstacles, which is at the borderline with image reconstructions with linear and nonlinear frameworks.²¹^,³⁷^–³⁹

Disclosures

The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.

Acknowledgments

The authors gratefully acknowledge the financial support of the Consejo Nacional de Ciencia y Tecnología (CONACYT), CONACYT-CICESE in La Paz (Mexico) proyecto interno No. 691114, within the program “Cátedras de Jóvenes Investigadores” and the Instituto de Ciencia y Tecnología del Distrito Federal (ICYT) in Mexico. This research was carried out with the administrative provision of the “Prosgrado en Ciencias e Ingeniería de la Computación” (PCIC) at the Universidad Nacional Autónoma de México (UNAM) and the support of the Seismic Laboratory for Imaging and Modeling (SLIM) at the University of British Columbia in Canada. The author CGS acknowledges the partial support from the programa PASPA-DGAPA at UNAM, and the BMO Institute at the Lübeck University for their host and support in sabbatical leave.

References

1.

M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 041101 (2006). https://doi.org/10.1063/1.2195024 RSINAK 0034-6748 Google Scholar

2.

M. K. A. Singh et al., “Photoacoustic reflection artifact reduction using photoacoustic-guided focused ultrasound: comparison between plane-wave and element-by-element synthetic backpropagation approach,” Biomed. Optics Express, 8 2245 –2260 (2017). https://doi.org/10.1364/BOE.8.002245 BOEICL 2156-7085 Google Scholar

3.

H.-M. Schwab, M. F. Beckmann and G. Schmitz, “Photoacoustic clutter reduction by inversion of a linear scatter model using plane wave ultrasound measurements,” Biomed. Opt. Express, 7 1468 –1478 (2016). https://doi.org/10.1364/BOE.7.001468 BOEICL 2156-7085 Google Scholar

4.

B. T. Cox and B. E. Treeby, “Artifact trapping during time reversal photoacoustic imaging for acoustically heterogeneous media,” IEEE Trans. Med. Imaging, 29 387 –396 (2010). https://doi.org/10.1109/TMI.2009.2032358 ITMID4 0278-0062 Google Scholar

5.

X. L. Dean-Ben et al., “Statistical approach for optoacoustic image reconstruction in the presence of strong acoustic heterogeneities,” IEEE Trans. Med. Imaging, 30 401 –408 (2011). https://doi.org/10.1109/TMI.2010.2081683 ITMID4 0278-0062 Google Scholar

6.

J. Jithin et al., “Speed-of-sound compensated photoacoustic tomography for accurate imaging,” Med. Phys., 39 7262 –7271 (2012). https://doi.org/10.1118/1.4764911 MPHYA6 0094-2405 Google Scholar

7.

H. Grün et al., “Photoacoustic tomography of heterogeneous media using a model-based time reversal method,” Proc. SPIE, 6856 685620 (2008). https://doi.org/10.1117/12.763111 PSISDG 0277-786X Google Scholar

8.

P. L. Rivière, J. Zhang and M. Anastasio, “Image reconstruction in optoacoustic tomography for dispersive acoustic media,” Opt. Lett., 31 (6), 781 –783 (2006). https://doi.org/10.1364/OL.31.000781 OPLEDP 0146-9592 Google Scholar

9.

H. Ammari, “Photoacoustic imaging for attenuating acoustic media,” Lect. Notes Math., 2035 57 –84 (2012). https://doi.org/10.1007/978-3-642-22990-9 LNMAA2 0075-8434 Google Scholar

10.

C. Xu et al., “Investigation of a diffuse optical measurements-assisted quantitative photoacoustic tomographic method in reflection geometry,” J. Biomed. Opt., 17 061213 (2012). https://doi.org/10.1117/1.JBO.17.6.061213 JBOPFO 1083-3668 Google Scholar

11.

K. Maslov and L. V. Wang, “Photoacoustic imaging of biological tissue with intensity-modulated continuous-wave laser,” J. Biomed. Opt., 13 024006 (2008). https://doi.org/10.1117/1.2904965 JBOPFO 1083-3668 Google Scholar

12.

M. Eghtedari et al., “High sensitivity of in vivo detection of gold nanorods using a laser optoacoustic imaging system,” Nano Lett., 7 1914 –1918 (2007). https://doi.org/10.1021/nl070557d NALEFD 1530-6984 Google Scholar

13.

Z. Zha et al., “Biocompatible polypyrrole nanoparticles as a novel organic photoacoustic contrast agent for deep tissue imaging,” Nanoscale, 5 4462 –4467 (2013). https://doi.org/10.1039/c3nr00627a NANOHL 2040-3364 Google Scholar

14.

V. M. Moock et al., “Frequency analysis for an extended photoacoustic transport model,” Opt. Lett., 40 4030 –4033 (2015). https://doi.org/10.1364/OL.40.004030 OPLEDP 0146-9592 Google Scholar

15.

A. I. Nachman, J. F. Smith and R. C. Waag, “An equation for acoustic propagation in inhomogeneous media with relaxation losses,” J. Acoust. Soc. Am., 88 1584 –1595 (1990). https://doi.org/10.1121/1.400317 JASMAN 0001-4966 Google Scholar

16.

P. L. Rivière, J. Zhang and M. Anastasio, “Image reconstruction in optoacoustic tomography accounting for frequency-dependent attenuation,” in IEEE Nuclear Science Symp. Conf. Record, 1841 –1845 (2005). https://doi.org/10.1109/NSSMIC.2005.1596689 Google Scholar

17.

V. Barnard, “The inversion of the Boltzmann transport equation with medical optical tomography applications,” (2006). Google Scholar

18.

V. M. Moock et al., “Photoacoustic tomography with diffusion approximation,” in Proc. of the WCE, 1 –4 (2013). Google Scholar

19.

S. Arridge, “Optical tomography in medical imaging,” Inverse Prob., 15 R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

20.

A. Ishimaru, “Diffusion of a pulse in densely distributed scatterers,” J. Opt. Soc. Am., 68 1045 –1050 (1978). https://doi.org/10.1364/JOSA.68.001045 JOSAAH 0030-3941 Google Scholar

21.

E. Gutiérrez-Reyes et al., “Heat transport approximation in the thermoelastic expansion. How necessary is this simplification in the laser induced photothermal phenomenon?,” J. Phys. Commun., 1 1 –16 (2018). Google Scholar

22.

G. Ku and L. V. Wang, “Deeply penetrating photoacoustic tomography in biological tissues enhanced with an optical contrast agent,” Opt. Lett., 30 507 –509 (2005). https://doi.org/10.1364/OL.30.000507 OPLEDP 0146-9592 Google Scholar

23.

B. Reyes-Ramírez, C. García-Segundo and A. García-Valenzuela, “On the spectral response of thick piezoelectric capacitive sensors for arrays in imagenology applications,” Proc. SPIE, 9040 90401J (2014). https://doi.org/10.1117/12.2043789 PSISDG 0277-786X Google Scholar

24.

E. Gutiérrrez-Reyes et al., “Analysis of the transfer function for layered piezoelectric ultrasonic sensors,” AIP Adv., 7 065106 (2017). https://doi.org/10.1063/1.4985562 AAIDBI 2158-3226 Google Scholar

25.

A. García-Valenzuela, R. G. Barrera and E. Gutiérrez-Reyes, “Rigorous theoretical framework for particle sizing in turbid colloids using light refraction,” Opt. Express, 16 19741 –19756 (2008). https://doi.org/10.1364/OE.16.019741 OPEXFF 1094-4087 Google Scholar

26.

L. D. Landau and E. M. Lifshitz, Theory of Elasticity, 7 2nd ed.Pergamon Press, Oxford (1970). Google Scholar

27.

A. D. Poularikas, The Handbook of Formulas and Tables for Signal Processing, CRC Press LLC, Boca Raton (1999). Google Scholar

28.

E. L. Madsen et al., “Anthropomorphic breast phantoms for testing elastography systems,” Ultrasound Med. Biol., 32 857 –874 (2006). https://doi.org/10.1016/j.ultrasmedbio.2006.02.1428 USMBA3 0301-5629 Google Scholar

29.

R. Ma et al., “Non-invasive whole-body imaging of adult zebrafish with optoacoustic tomography,” Phys. Med. Biol., 57 7227 –7237 (2012). https://doi.org/10.1088/0031-9155/57/22/7227 PHMBA7 0031-9155 Google Scholar

30.

M. Heijblom et al., “Visualizing breast cancer using the Twente photoacoustic mammoscope: what do we learn from twelve new patient measurements?,” Opt. Express, 20 11582 –11597 (2012). https://doi.org/10.1364/OE.20.011582 OPEXFF 1094-4087 Google Scholar

31.

G. T. Herman, Fundamentals of Computerized Tomography: Image Reconstruction from Projections, Springer, London (2009). Google Scholar

32.

G. T. Herman and D. Odhner, “Performance evaluation of an iterative image reconstruction algorithm for positron emission tomography,” IEEE Trans. Med. Imaging, 10 336 –346 (1991). https://doi.org/10.1109/42.97583 ITMID4 0278-0062 Google Scholar

33.

S. Telenkov and A. Mandelis, “Signal-to-noise analysis of biomedical photoacoustic measurements in time and frequency domains,” Rev. Sci. Instrum., 81 124901 (2010). https://doi.org/10.1063/1.3505113 RSINAK 0034-6748 Google Scholar

34.

D. J. Schroeder, Astronomical Optics, Academic Press, San Diego (1999). Google Scholar

35.

J. Pinkerton, “The absorption of ultrasonic waves in liquids and its relation to molecular constitution,” Proc. Phys. Soc. B, 62 129 –141 (1948). https://doi.org/10.1088/0370-1301/62/2/307 PPSBAP 0370-1301 Google Scholar

36.

L. D. Landau and E. M. Lifshitz, Fluid Mechanics, 2nd ed.Butterword and Heinemann, Oxford (1987). Google Scholar

37.

C. Buj et al., “Noncontact holographic detection for photoacoustic tomography,” J. Biomed. Opt., 22 106007 (2017). https://doi.org/10.1117/1.JBO.22.10.106007 JBOPFO 1083-3668 Google Scholar

38.

S. Marburg and B. Nolte, Computational Acoustics of Noise Propagation in Fluids, Springer-Verlag, Berlin (2008). Google Scholar

39.

F. Filbir et al., in Modeling, analysis, and approximation theory toward applications in tomography and inverse problems, (2018). Google Scholar

Biography

Verena Margitta Moock received her PhD in computer science from the Universidad Nacicional Autónoma de México (UNAM) in 2015 and her MSc and BSc degrees in applied mathematics from the Technische Universität Darmstadt. She is involved with research on photoacoustics at the Instituto de Ciencias Aplicadas y Tecnología (ICAT) at UNAM. Since 2017, she forms part of the National System of Researchers (SNI) of the Consejo Nacional de Ciencia y Tecnología (CONACYT) of Mexico and is a lecturer at Colegio Alemán Alexander von Humboldt A.C. in Mexico City. Between 2016 and 2017 she went to the Seismic Laboratory for Imaging and Modeling of the University of British Columbia in Vancouver, Canada, for a postdoctoral visit.

Edahí Antonio Gutiérrez Reyes received his DSc degree from the Universidad Nacional Autónoma de México (UNAM) in 2013 and his BSc degree in physics from UNAM in 2005. He worked 1 year in a postdoctoral fellowship at Benemérita Universidad Autónoma de Puebla in 2014 and in a postdoctoral fellowship at CCADET UNAM in 2015. Now he works at centro de investigación científica y de educación superior de Ensenada, Baja California, unidad la Paz CICESE. He works on seismic hazard, electrodynamics, photoacoustics, elasticity, and wave propagation in heterogeneous systems.

Crescencio Garcia-Segundo received his PhD in optical sciences from the Instituto Nacional de Astrofísica, Óptica y Electrónica, Mexico, in 1996. His research focuses on photoacoustics, photothermal imaging, and biophotonics in random media. He is a tenure researcher at ICAT and a member of the National System of Researchers of Mexico. He has collaborations with the Institute of Biomedical Optics at the University of Lübeck, at the Wuhan Institute of Physics and the Imperial College London.

Citation Download Citation

Verena Margitta Moock, Edahí Antonio Gutiérrez-Reyes, and Crescencio García-Segundo "Image reconstruction with the Heaviside equation in photoacoustic tomography accounting for dispersive acoustic media," Journal of Biomedical Optics 23(7), 076010 (19 July 2018). https://doi.org/10.1117/1.JBO.23.7.076010

Received: 12 December 2017; Accepted: 21 June 2018; Published: 19 July 2018

Access the abstract

JOURNAL ARTICLE
12 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 1 scholarly publication.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Acoustics

Image restoration

Signal attenuation

Inverse problems

Photoacoustic tomography

Signal to noise ratio

Wave propagation

1.

Introduction

2.

Inverse Problems in Waveform Tomography