Modeling photoacoustic imaging with a scanning focused detector using Monte Carlo simulation of energy deposition

Guenther Paltauf; Paul R. Torke; Robert Nuster

doi:10.1117/1.JBO.23.12.121607

24 September 2018 Modeling photoacoustic imaging with a scanning focused detector using Monte Carlo simulation of energy deposition

Guenther Paltauf, Paul R. Torke, Robert Nuster

Author Affiliations +

Journal of Biomedical Optics, Vol. 23, Issue 12, 121607 (September 2018). https://doi.org/10.1117/1.JBO.23.12.121607

Abstract

Photoacoustic imaging using a focused, scanning detector in combination with a pulsed light source is a common technique to visualize light-absorbing structures in biological tissue. In the acoustic resolution mode, where the imaging resolution is given by the properties of the transducer, there are various challenges related to the choice of sensors and the optimization of the illumination. These are addressed by linking a Monte Carlo simulation of energy deposition to a time-domain model of acoustic propagation and detection. In this model, the spatial and electrical impulse responses of the focused transducer are combined with a model of acoustic attenuation in a single response matrix, which is used to calculate detector signals from a volumetric distribution of absorbed energy density. Using the radial symmetry of the detector, the calculation yields a single signal in less than a second on a standard personal computer. Various simulation results are shown, comparing different illumination geometries and demonstrating spectral imaging. Finally, simulation results and experimental images of an optically characterized phantom are compared, validating the accuracy of the model. The proposed method will facilitate the design of photoacoustic imaging devices and will be used as an accurate forward model for iterative reconstruction techniques.

1. Introduction

Photoacoustic microscopy (PAM) is a well-established method for generating images of biological tissue containing structures with optical absorption contrast. In the acoustic resolution mode, where the lateral resolution is given by the focusing properties of an ultrasound transducer that scans across the tissue surface, an imaging depth of several millimeters has been observed.¹ High repetition-rate pulsed laser systems and fast scanning devices allow short image acquisition times, despite the necessity for moving the transducer point by point over the tissue surface. In this study, we use the term PAM to describe a method where a scanning, focused detector receives signals from sources along its line of sight, and the acquired signals (A-scans) are combined to a two- or three-dimensional (2-D or 3-D) image. This includes scenarios for which sometimes also the terms “mesoscopy” or “macroscopy” are used to indicate that the observed objects may actually not be in the micrometer-range but can reach sizes on the order of millimeters.²^,³ By contrast, optical resolution PAM can image submicrometer structures within a depth of $\sim 1 mm$ , only relying on the optical focusing of the illuminating lens.⁴

Accurate modeling is essential in the design of a PAM system for a specific application. By employing a realistic imaging model, factors such as illumination geometry and various transducer properties can be optimized. “Realistic” means that as many factors as possible, which play a role in the formation of an image, are taken into account. Specifically, the problem that acoustic and optical paths have to be separated in a PAM device can be addressed in a model that includes both, the optical illumination and the acoustic detection processes. The optical part has to provide accurate predictions of the energy deposition map for spatially heterogeneous distributions of optical properties in the object and for various illumination patterns on the surface of the object. The acoustical part should model the generation, propagation, and detection of photoacoustically generated acoustic waves. This includes acoustic damping effects and the integration of the waves over the finite area of an ultrasound detector.

Concerning the optical illumination in PAM, the method of choice for modeling light transport in biological tissue is the Monte Carlo (MC) simulation.⁵ There are several examples of how MC simulations have been used to optimize the light distribution for photoacoustic imaging, for instance, for obtaining a maximum signal from a source at a given depth⁶ or for optimizing dark field illumination in PAM.⁷ For a complete modeling of the photoacoustic imaging process, results of the MC simulation can be coupled with the acoustic generation and propagation. Jacques⁸ used MC-generated energy distribution maps to simulate signals of point-like detectors. Ermilov et al.⁹ used MC results to model signals received by an array of flat, finite-sized detector elements.

A well-established method for simulating the signal generation in photoacoustic imaging is the k-wave toolbox for MATLAB.¹⁰ It can model acoustic wave generation from an arbitrary distribution of absorbed energy density and is able to take into account heterogeneous distributions of sound speed. MC simulations combined with k-wave were used by Kaplan et al.¹¹ for quantitative photoacoustic imaging or by Heijblom et al.¹² to improve interpretation of clinical photoacoustic images. Additionally, k-wave has been used to model PAM with optical resolution, recording acoustic waves by a point-like transducer that scans across an object.¹³^,¹⁴ The pseudospectral method computes the evolution of the sound field in the entire source volume at a series of time steps. Multiple detector positions can be defined within the computational grid, at which the temporal pressure signals are recorded. This feature also allows simulating an extended detector surface by summing point detector signals; however, 3-D simulations with the k-wave toolbox are quite time consuming. For instance, a single time step in a computational grid with a size of $256^{3} voxels$ requires about one second on the CPU of a multicore desktop computer.¹⁰ Moreover, in the simulation of PAM with acoustic resolution (AR-PAM) the volume containing the photoacoustic sources and the focusing detector are separated by a distance equal to the acoustic focal length. If both have to be modeled with sufficiently fine computational grids to resolve acoustic frequencies in the MHz range, the computational domain in the k-wave simulation would become impractically large.

An alternative is time-domain modeling of acoustic propagation and detection, using the spatial impulse response (SIR), which describes the temporal spreading of an acoustic wave from a point source over the detector surface. Various implementations of the SIR have been used to improve reconstructions in photoacoustic tomographic imaging.¹⁵^,¹⁶ In PAM, Araque et al.¹⁷ demonstrated an improved depth of field, if a model-based reconstruction applied to line-scanning detector data used the SIR of a focused detector. The SIR was used by Karmakar et al.¹⁸ to simulate the response of a focused detector to a distribution of point-like acoustic sources. The concept of the SIR is also widely applied in simulations of pulse-echo ultrasound imaging, where it is implemented in the FIELD II software.¹⁹^,²⁰ The amount of data contained in the SIR can be huge. For an arbitrary shape of the detection surface, the signal from every point within the volume containing the acoustic sources can be different. If precomputed SIR data are stored in a system matrix, the required storage can, therefore, become very large.

As the SIR concept is well suited to model acoustic wave generation and detection in AR-PAM, we developed a simulation method that avoids large storage requirements by taking advantage of the radial symmetry of focusing sensors. Furthermore, the computation time for a single acoustic signal on the CPU of a standard PC is shown to be less than a second even if it involves integration over $256^{3}$ grid points of the MC-simulated energy density map. The model also includes a power law acoustic attenuation and the finite bandwidth of the detector.

In the following sections, we start with a short description of the required specifications for the light transport model. Then, we proceed with an accurate explanation of the acoustic part of the model, including the calculation of the sensor SIR, the incorporation of acoustic damping and detector bandwidth, and the scanning of the detector. Simulation results are then presented, where we compare dark field and bright field illumination techniques in combination with focusing sensors. A second set of simulations demonstrates spectrally resolved PAM of blood vessels with different oxygenation levels. Finally, we show a comparison between simulation and experiment, using a ring-shaped sensor that provides a large depth of field.

2. Methods

2.1.

Monte Carlo Simulation of Light Transport

The simulation of sound waves generated by short laser pulses via the photoacoustic effect in biological tissue requires the knowledge of a distribution of absorbed energy density $W$ , which is given by the fluence $Ψ$ and the absorption coefficient $μ_{a}$ . From the energy density distribution, the initial pressure $p_{0}$ is derived by multiplication with the Gruneisen parameter $Γ$

Eq. (1)

p_{0} (r) = Γ (r) W (r) = Γ (r) μ_{a} (r) Ψ (r) .

An MC simulation for optically heterogeneous media that models light and deposited energy distributions in biological tissue has been developed by Jacques et al.⁸^,²¹ The 3-D simulation assigns values of $μ_{a}$ , the scattering coefficient $μ_{s}$ , and of the scattering anisotropy $g$ to each voxel in a volumetric Cartesian grid. The wavelength-dependent values of the optical properties are specific for different types of tissues.²² As the MC simulation is a standard method, it will not be described in detail here. Only the requirements for the further simulations of photoacoustic signals will be specified. In particular, the distributions of $W$ and $μ_{a}$ should be stored on a sufficiently fine grid, for instance on $256^{3}$ grid points. The grid spacing determines the achievable acoustic frequency as described in the next section. In our simulations, we assumed $Γ$ constant. If this is not the case, its distribution should be stored as well. Furthermore, the laser radiation should be delivered in a beam of finite diameter. In this study, we used a matched boundary, where the refractive index of the tissue and the adjacent medium are identical. This is a good approximation for the case where an acoustically matching material, such as water or ultrasound gel is applied to the tissue surface for acoustic coupling to the detector.

2.2.

Photoacoustic Sound Generation

First, the continuous model for photoacoustic sound generation from a distribution of initial pressure and its reception by a sensor with finite area is described, followed by the discrete implementation of the model. Given a distribution of absorbed energy density $W (r)$ , the time-dependent pressure signal $p (r_{0}, t)$ at a detection point $r_{0}$ is given by the following integral:¹⁵

Eq. (2)

p (r_{0}, t) = \frac{β}{4 π C_{p}} \frac{\partial}{\partial t} \int_{V} W (r) \frac{1}{| r_{0} - r |} δ (t - \frac{| r_{0} - r |}{c_{s}}) d^{3} r,

where

c_{s}

is the speed of sound,

β

is the thermal expansion coefficient,

C_{p}

is the specific heat capacity at constant pressure, and

V

is the volume that contains the

W

-distribution. Focused detection requires a sensor with finite area, which obtains its focusing property from a specific shape. Most common is a spherically focused sensor, where the center of curvature determines the focus position. Also possible are annular sensors, which focus along the symmetry axis of the ring.³^,²³^,²⁴ Similar focusing properties can be achieved with a conical sensor, which also preferentially receives signals from its symmetry axis.²⁵ Our analysis will be limited to all kinds of detectors with radial symmetry, as only these types of detectors exhibit a focus or a focal range suitable for scanning photoacoustic imaging. The signal received by a sensor with finite area is calculated by integrating the point detector signal from Eq. (2) over the sensor surface

Eq. (3)

s (r_{d}, t) = \int_{S_{0}} p (r_{0}, t) d S_{0} (r_{d}) = \frac{β}{4 π C_{p}} \frac{\partial}{\partial t} \int_{S_{0}} \int_{V} W (r) \frac{1}{| r_{0} - r |} δ (t - \frac{| r_{0} - r |}{c_{s}}) d^{3} r d S_{0} (r_{d}),

where

S_{0}

denotes the sensor surface and

r_{d}

is the position of the center of the sensor. Separating the source

W

from the sensor-dependent part of this integral yields

Eq. (4)

s (r_{d}, t) = \int_{V} W (r) h^{i d} (r_{d}, r, t) d^{3} r .

The function

h^{i d} (r_{d}, r, t)

describes the system response of the sensor to a photoacoustic point source

Eq. (5)

h^{i d} (r_{d}, r, t) = \frac{β}{4 π C_{p}} \frac{\partial}{\partial t} \int_{S_{0}} \frac{1}{| r_{0} - r |} δ (t - \frac{| r_{0} - r |}{c_{s}}) d S_{0} (r_{d}) .

The ideal system response

h^{i d} (r_{d}, r, t)

defined here depends only on the position of a point in the source volume relative to the sensor and assumes an ideal, lossless medium and a detector with uniform frequency response. To take into account the acoustic damping and the specific electrical response of the detector, a more general definition of the system response can be used, where

h

is now a convolution of several response functions

Eq. (6)

h (r_{d}, r, t) = g^{att} (r_{d}, r, t) * g^{e l} (t) * g^{p a} (t) * g^{s} (r_{d}, r, t),

where “

*

;” describes temporal convolution,

g^{e l}

is the acoustoelectrical impulse response (EIR) of the detector,

g^{a t t}

is the attenuation impulse response,

g^{s}

is the SIR, and

g^{p a}

describes the specific temporal signal from a photoacoustic point source. The impulse response defined in Eq. (5) can be written as

h^{i d} (r_{d}, r, t) = g^{p a} (t) * g^{s} (r_{d}, r, t)

, where

g^{s}

describes the spreading of a delta pulse over the sensor surface and

g^{p a}

is a time derivative operator. If the laser pulse has a finite duration, then

g^{p a}

is the time derivative of the temporal pulse shape. It is possible to derive analytical equations for the SIR in various geometries.²⁶ In the discrete implementation we use a numerical calculation, which can be easily adapted for a variety of transducer types, such as flat, annular, or spherically focused.

We proceed with the discrete model of photoacoustic sound generation and detection and start with the signal on a single position. Scanning the detector will be described later. Analogous to Eq. (4), the signal vector $s$ with components $s_{k} \equiv s (t_{k})$ , where $t_{k} = k Δ t$ is generated by an absorbed energy density distribution contained in a vector $W$ with components $W_{j} \equiv W (r_{j})$ , where $r_{j}$ defines a point on a 3-D Cartesian grid with spacings $Δ x$ , $Δ y$ , and $Δ z$ , and is given as

Eq. (7)

s_{k} = \sum_{j} h_{k j} W_{j} or s = H W .

Spatial and temporal increments depend on each other. According to the sampling theorem, the minimum acoustic wavelength generated with a grid spacing

Δ x

(assuming

Δ x = Δ y = Δ z

) is

λ_{\min} = 2 Δ x

, leading to a maximum acoustic frequency

f_{\max} = c_{s} / (2 Δ x)

. To avoid a loss of information, the temporal increment

Δ t

should not be larger than

Δ x / c_{s}

. The signal vector

s

has

K

components and the source distribution is defined on

J

discrete positions

r_{j}

relative to the sensor. The system response

H

with components

h_{k j}

is, therefore, a matrix of size

K \times J

. Also, the SIR,

G^{s}

, has the same size, whereas the discrete counterparts of

g^{a t t}

,

g^{e l}

, and

g^{p a}

are the square matrices

G^{a t t}

,

G^{e l}

, and

G^{p a}

of size

K \times K

.

For the discrete, numerical calculation of the SIR, the curved surface of the detector is first divided into surface elements located at $r_{0, i}$ . The spherical wave emitted by the point source hits different locations of the detection surface at different times. As the SIR is a function of time, the task is to find at each specific time the number and strength of spherical waves hitting the surface elements and summing up all these contributions. This corresponds to the integral over the delta function in Eq. (5). The components of $G^{s}$ , $g_{k j}^{s}$ are given as

Eq. (8)

g_{k j}^{s} = \sum_{i} w_{i} q_{k i j}, with q_{k i j} = {\begin{cases} \frac{1}{| r_{0, i} - r_{j} |} & if | \frac{| r_{0, i} - r_{j} |}{c_{s}} - t_{k} | \leq \frac{Δ t}{2} \\ 0 & otherwise \end{cases} .

The weight factors

w_{i}

take into account the size of surface elements at

r_{0, i}

. The next step is the definition of the photoacoustic signal

G^{p a}

. Commonly, a photoacoustic point source is modeled as a solid sphere, for which the analytical signal has a typical “

N

”-shape.¹⁵^,²⁷ Defining the sphere radius

a = Δ x

, where

Δ x

is the spatial increment of the

W

-distribution, gives

Eq. (9)

g_{k + m, k}^{p a} = {\begin{cases} - \frac{c_{s}}{2 a} m Δ t & if ? - \frac{a}{c_{s} Δ t} \leq m \leq \frac{a}{c_{s} Δ t} \\ 0 & otherwise \end{cases},

where

m

is an integer. The resulting square matrix has the structure of a Toeplitz matrix, with replications of the basic “

N

”-signals grouped in columns along its main diagonal. The temporal convolution with

G^{s}

is then achieved by a matrix multiplication. It is possible to define convolution matrices

G^{a t t}

and

G^{e l}

separately. For instance, in some cases the EIR of the sensor is known and can be implemented in a convolution matrix in a similar way as the photoacoustic signal matrix

G^{p a}

was obtained. Sometimes the detector transfer function is given or its center frequency and a fractional bandwidth are known. This suggests defining a suitable filter in frequency space that can be multiplied with the Fourier transform of the photoacoustic response function in the columns of

G^{p a}

. The result is a modified matrix

{\tilde{G}}_{p a}

that now also contains the EIR. It is important to ensure causality of the obtained filter. In other words, the simulated sensor should not start to respond before the actual sound wave has arrived, as it would be the case for a noncausal filter. If the amplitude spectrum is given, the related phase spectrum for a minimum phase causal filter can be found by applying the Hilbert transform to the logarithm of the amplitude spectrum.²⁸

In the same way as the detector response function, also the attenuation can be taken into account by defining a suitable filter and applying it to the photoacoustic response function in frequency space. The difference is that the attenuation is not a function of frequency alone but also depends on the propagation distance of the wave. The amplitude attenuation is usually expressed by a power law

Eq. (10)

g^{a t t} (f, | r_{0} - r |) = \exp (- α | r_{0} - r |) with α = α_{0} f^{b},

f

is the frequency,

α

is the acoustic absorption coefficient, and

α_{0}

is a constant. This kind of filter does not take into account the dispersion of the wave and therefore leads to a noncausal impulse response. For biological tissue

b \approx 1

, for which case a minimum phase causal filter that takes into account dispersion was proposed by Gurumurthy²⁸

Eq. (11)

g^{a t t} (f, | r_{0} - r |) = \exp (- α_{0} | r_{0} - r | f) \cdot \exp [i \frac{2 f}{π} α_{0} | r_{0} - r | \ln | \frac{f}{f_{0}} |] .

A value of

f_{0} = 77 MHz

was suggested in Ref. 28 [actually the authors introduce a minimum phase delay factor

τ_{m} = \ln (2 π f_{0}) = 20

]. Instead of this relation, any other causal impulse response can be used, also for

b \neq 1

.

For a finite size detector, the distance from a source point to different parts of the detector surface is generally not constant; it is, therefore, not possible to define a unique filter function. This can be solved by an approximation, where an average distance is used for defining the attenuation. Another possibility is nonstationary filtering, where the filter function changes as a function of time.²⁹ Its implementation is straightforward as every column of $G^{p a}$ corresponds to a certain time of flight $t = | r_{0} - r | / c_{s}$ . Substituting $| r_{0} - r |$ by $c_{s}$ $t$ in Eq. (11) gives a slightly different filter for each column in the Fourier transformed matrix $G^{p a}$ and consequently different kernels in each of its columns after inverse Fourier transform. This results in a final, modified photoacoustic response matrix ${\tilde{G}}_{p a}$ that leads to a nonstationary convolution when it is multiplied with $G^{s}$ . An alternative to the structure of this convolution matrix is “nonstationary combination,” as it was suggested by Margrave²⁹ and used by Treeby for building a convolution matrix that compensates attenuation in a photoacoustic imaging reconstruction algorithm.³⁰ In this case, the filter kernel for a given time and thus propagation distance is written in the row of the filter matrix, therefore describing a filter that changes as a function of the filter output time and can take into account abrupt changes of the filter characteristics.

Due to the radial symmetry of the sensor, it is convenient to use cylindrical coordinates for the calculation of the SIR (see Fig. 1), which, therefore, becomes a function $g^{s} [ρ (r_{d}, r), ζ (r_{d}, r), t]$ . As the MC simulation usually gives a $W$ -distribution in Cartesian coordinates, there are two possibilities to implement the summation in Eq. (7). Either the summation goes over the Cartesian grid, after calculating the position $(ρ, ζ)$ of each grid point relative to the detector. This results in a relatively slow summation over the source volume. The second option is an interpolation of the $W$ -distribution into cylindrical coordinates $(ρ, ϕ, ζ)$ , followed by an integration over $ϕ$ . While the detector is scanned, the origin of the cylindrical coordinate system, $r_{d}$ , moves relative to $W$ . The summation in Eq. (7) is then performed over the two remaining coordinates $ρ$ and $ζ$ , respectively. As indicated above, this is done by multiplying the system response matrix with the vector containing the source distribution. Overall, in our implementation in MATLAB (R2015b) on a standard PC (quad-core 3.2GHz CPU), the second option resulted in about 50 times faster computation compared with the first one. With this acceleration, the calculation of a single A-scan signal took about 0.5 s.

Fig. 1

Arrangement of detector and absorbed energy density distribution as it is used to calculate the detector SIR. The energy density is given in Cartesian coordinates $(x, y, z)$ .The detector is defined in a cylindrical coordinate system $(ρ, ζ)$ with its origin at $r_{d}$ . Spherical acoustic waves generated at points $r$ are received by points $r_{0}$ on the detector surface.

2.3.

Scanning the Detector

The next step is the simulation of B-scan images. With the precomputed SIR of the sensor, the time required to calculate a 2-D B-scan consisting of 100 A-scans from a volume containing $256^{3}$ entries of absorbed energy density is about one minute. However, this assumes that the energy density distribution is also precomputed and does not need to be repeatedly calculated during the simulation.

This assumption is used in the first described method: the illumination of the sample remains static while the focused sensor moves relative to the constant $W$ -distribution. An imaging scenario, where such an arrangement of illumination and detection would be applicable, is PAM in the transmission mode. An object is illuminated with laser pulses from one side and a focused sensor scans the surface of the object from the other side, recording acoustic waves that have propagated in the direction of the incoming light.²

In many imaging scenarios, however, illumination and collection of acoustic signals have to take place on the same side of the object (“reflection mode”). In this case, the illumination of the object generally moves with the sensor and does not stay constant. This leads to the second method, where at each new sensor position a new MC simulation has to be run. An example is acoustic resolution photoacoustic microscopy (AR-PAM) with dark field illumination, where special optics creates a ring-shaped illumination pattern around the ultrasound sensor.¹ In addition to bypassing the optically opaque piezoelectric sensor, this kind of illumination also diminishes the strong signal generated at the surface. Another example is AR-PAM with bright field illumination through a hole in the center of a focusing detector.³¹ This method accepts some disturbance due to the strong surface signal but benefits from a much simpler design of the illumination optics. To study the differences between these illumination modes, the influence of the light distribution moving with the sensor cannot be neglected. If the MC simulation is the computational bottleneck, calculation of a new energy density map at each scanning position will drastically increase the overall simulation time. This approach may be feasible when performing the MC simulation on a graphics processing unit, where a significant acceleration compared with a standard processor is reported.³²^,³³

Finally, a method is presented that avoids calculation of multiple MC simulations but still allows taking into account the changes in the energy density distribution when the illumination moves with the detector. It gives a reasonable estimation if the embedded structures have relatively low absorption. Such a situation may arise in the near infrared, where blood vessels still have good contrast in PA images, but the blood absorption is much smaller than in the visible part of the spectrum. It can be assumed that within the weakly absorbing structure the fluence remains at its undisturbed value $Ψ_{0}$ , calculated with the background $μ_{a}$ and $μ_{s}$ values. In a first step, this fluence is calculated for a static illumination of the phantom. To simulate a B-scan, the sensor is positioned at its designated location relative to the illumination beam and the $μ_{a}$ -distribution of the embedded absorbing structures, such as blood vessels, is shifted across the fluence distribution, calculating at each position $W_{approx} = μ_{a}$ $Ψ_{0}$ . Figure 2(b) shows a section through an energy density distribution of a phantom containing several absorbing objects. The values of absorption and scattering coefficient for these objects are shown in Fig. 2(a). The background properties in this phantom are from breast tissue, covered by a skin (epidermis) layer containing melanin.²² The phantom was illuminated with a flat top beam with 0.6-cm diameter. Figures 2(c) and 2(d) show profiles through the energy density distribution, comparing at different positions in the phantom the approximation outlined above ( $W_{approx}$ ) with the full MC simulation ( $W_{exact}$ ), which correctly models the drop of the fluence inside an absorbing object. The best coincidence is achieved for the spherical object with the lowest absorption. If the absorption becomes higher, such as in the simulated blood vessels, the approximation tends to overestimate the energy density. This is expected as the approximation assumes the undisturbed fluence inside the absorbing structure. Nevertheless, it is reasonably accurate and allows a comparison between different illumination geometries, assuming that the relative error due to the overestimated fluence remains similar.

Fig. 2

(a) Outline of the numerical phantom containing several cylindrical blood vessels with 0.05- and 0.1-cm diameter, one “tumor” with 0.2-cm diameter, and an epidermal layer containing melanin ( $120 - μ m$ thickness). Optical properties at 750-nm wavelength, with an anisotropy factor $g = 0.9$ for the whole phantom. (b) Logarithmic display of the energy density in the center plane of the phantom after illumination with a 0.6-cm diameter beam from below. (c) Profile through the exact energy density distribution $W_{exact}$ along line A, compared with the approximation $W_{approx} = μ_{a} Ψ_{0}$ . (d) Profiles along line B.

2.4.

Experiment

An experiment was designed to compare the output of the model with the real-world data. In this experiment, it was important to build a phantom with well-defined optical properties, which could be used as an input for the MC simulation. The outline of the phantom and the sensor is shown in Fig. 3. The phantom consisted of a block of agar (2% agar in water) containing 4% SMOFlipid (20% fat content, similar to Intralipid). A small amount of indocyanine green (ICG) at a concentration of about $0.5 μ g / mL$ was added to the agar to achieve a background absorption coefficient of $0.15 {cm}^{- 1}$ at 750-nm wavelength. The reduced scattering coefficient of the agar was $μ_{s}^{'} = 6 {cm}^{- 1}$ , measured with oblique incidence video reflectometry.³⁴ The absorbing structures within the agar phantom were silicon tubes with an inner diameter of 0.1 cm, filled with a solution of ICG with a concentration of $90 μ g / mL$ and $μ_{a} = 20 {cm}^{- 1}$ , arranged at depths of 0.5 and 0.75 cm, respectively. The phantom was scanned with a detector made of the piezoelectric polymer PVDF. The sensor had a thickness of $110 μ m$ and had the shape of a ring with an inner radius of 1.0 cm and an outer radius of 1.18 cm. The annular piece of PVDF was glued onto a conical substrate, which was inclined toward the ring axis by an angle of 25 deg. Such a sensor has the property to focus onto a line, providing a large depth of field. It combines properties of a large area conical sensor²⁵ and a flat ring,²⁴^,²³ and is currently under investigation as an element of an array of flat and inclined ring sensors.³⁵ The spectral response of the PVDF sensor was modeled using a method originally proposed by Cox and Beard³⁶ to simulate the frequency-dependent directivity of an optical polymer film used as Fabry–Perot interferometric ultrasound sensor. It calculates the spectral response to a plane wave incident on a three-layer system consisting of water, the piezoelectric PVDF, and a matched backing material with the same properties as PVDF. The output from an optical parametric oscillator, pumped by a frequency-doubled Nd:YAG laser with 10-ns pulse duration, was guided to the sensor through a 1-mm core diameter optical fiber. The end face of the fiber was imaged with a lens through the center hole of the ring sensor onto the surface of the phantom, generating a homogeneously illuminated area with 0.6-cm diameter. In the simulation, the approximation was used where first the fluence distribution with fixed illumination of a phantom without absorbers was calculated and then the absorbing structures were scanned across the sensor. Although the resulting energy density, given as the product of absorption coefficient and undisturbed fluence, certainly overestimated the actual values, the temporal signals and the relative amplitudes of waves from different depths could be quite accurately modeled in this way. As the acoustic waves in this experiment propagated primarily in water, the acoustic damping properties of water were used in the simulation, with an attenuation coefficient that depends on frequency squared.³⁷

Fig. 3

Layout of the phantom and the sensor for the experiment. The phantom is a block of light-scattering agar containing two plastic tubes filled with an ICG solution, illuminated from below with laser pulses. Acoustic waves are recorded with a ring-shaped, $110 - μ m$ thick PVDF film detector mounted on a conical plastic block.

For the ring sensor, there is no linear relationship between the depth of the photoacoustic source and the time of flight of the wave. To reconstruct depth profiles, we multiplied the transposed impulse response matrix $H^{T}$ , calculated only for points with $ρ = 0$ , with the temporal signal vector $s$ . The effect of applying $H^{T}$ is in this case a transformation from temporal to spatial profiles along the detector axis. For this reconstruction, $H$ was derived from the pure SIR of the ideal sensor, without the finite sensor bandwidth and without acoustic damping.

3. Results

3.1.

Illumination Modes

In these simulations, the influence of the illumination geometry on B-scan images of the phantom shown in Fig. 2(a), recorded with a spherically focusing detector was investigated. The detector was modeled with a center frequency of 5 MHz, 100% bandwidth, a focal length of 2 cm, and a numerical aperture of 0.5. One purpose of this simulation was the investigation of deep imaging in the near infrared. The focus of the sensor was, therefore, located at 1-cm depth below the tissue surface. The acoustic damping of tissue calculated with Eq. (11) was used for the system response. Figure 4 shows how the EIR of this sensor, together with the acoustic attenuation, change the photoacoustic signal $g^{p a}$ from the small, spherical source. The figure compares $g^{p a}$ to the convolution $g^{a t t} (r_{d}, r, t) * g^{e l} (t) * g^{p a} (t)$ for two different values of $| r - r_{0} |$ , 1 and 3 cm, denoting the distance between the photoacoustic point source and a point on the sensor surface. These two signals are two columns of the modified convolution matrix ${\tilde{G}}_{p a}$ defined above, which performs the nonstationary filtering.

Fig. 4

Comparison of the photoacoustic signal $g^{p a}$ defined in Eq. (9) from a spherical source with $117 - μ m$ radius with the convolution $g^{a t t} * g^{e l} * g^{p a}$ for different values of the propagation distance $| r - r_{0} |$ : (a) 1 cm and (b) 3 cm.

The layout of the energy density distribution and the sensor is shown in Figs. 5(a) and 5(d), at a position where the detector is centered relative to the phantom. In this simulation, the spatial increment in the numerical phantom with a size of $3 \times 3 \times 3 c m^{3}$ was $Δ x = 117 μ m$ . This leads to a maximum generated acoustic frequency of $f_{max} = 6.4 MHz$ . The temporal increment was $Δ t = 40 ns$ , leading to a sampling rate two times the minimum required for accurately resolving features with $f_{max}$ . The size of the impulse response matrix was given by the number of temporal samples, $K = 976$ , multiplied by the number of source positions $(ρ, ζ)$ in the cylindrical coordinate system of the sensor, $J = 182 \times 257 = 46,774$ . The resulting double precision matrix occupied about 350 MB of storage.

Fig. 5

Comparison between dark field and center illumination. (a), (b), and (c): Dark field illumination; (d), (e), and (f): center illumination. (a) and (d) Outline of energy density distribution and focusing sensor. (b) and (e) B-scan images, (c) and (f) profiles along the dotted lines in (b) and (e).

For the dark field illumination, a ring with an inner radius of 0.6 cm and an outer radius of 1 cm was illuminated with a fluence corresponding to the maximum permissible exposure (MPE) at the selected wavelength of 750 nm, $H_{0} = 25.2 mJ / c m^{2}$ . The incident angle of the annular beam was 45 deg. In this particular geometry, the incident beam would just pass outside the acoustic detector and would not create any strong surface signals within the detector’s receiving aperture. For the center illumination, the fluence distribution was calculated by assuming a beam diameter of 0.3 cm at the tissue surface, again with the same MPE value. This beam diameter simulates illumination through a center hole in the sensor (1-mm radius) with an optical fiber having an NA of 0.4. As the illumination had to move with the sensor, the approximation was used in which the absorbing structures move across the fluence distribution and at each scanning position the energy density is calculated as $W_{approx} (r) = μ_{a} (r) Ψ_{0} (r)$ . The epidermis layer was present at all scanning positions and was, therefore, included in the calculation of $Ψ_{0}$ . Figures 5(b) and 5(e) show the B-scans for each of the two cases. These images display the envelope of the pressure signals, calculated with the help of the Hilbert transform. In both cases, the signal from the surface is rather high and is therefore cropped, providing better visibility of the deep structures. Although the total pulse energy is higher for the axicon illumination, the relative amplitude of the surface signal is similar (0.4) in both cases, as seen in the profiles displayed in Figs. 5(c) and 5(f). The main difference is the better visibility of the deep structures for the dark field illumination, where the ratio of signal amplitudes from buried sources to the surface signal is much higher. The dark field illumination, on the other hand, creates strong acoustic sources, which lie outside the focal volume of the transducer but are still detected. This is a kind of “clutter,” which manifests itself as a hazy background in Fig. 5(b).

3.2.

Spectrally Resolved Imaging

This simulation used a forward signal acquisition geometry, where the laser pulses irradiated a tissue slab from one side and the detector scanned the object from the opposite side. The illumination with a 0.6-cm diameter beam stayed at a fixed position relative to the phantom and the spherically focused detector moved relative to the obtained energy density distribution. The purpose of this simulation was to analyze the wavelength-dependent photoacoustic signal amplitudes of vessels containing blood with different oxygenation levels. For this analysis, we chose the two wavelengths $λ_{1} = 797 nm$ , where the absorption coefficients of oxygenated and deoxygenated hemoglobin are equal, and $λ_{2} = 750 nm$ , where deoxygenated Hb absorbs more strongly. Figure 6(a) shows the absorbed energy density in the phantom and indicates the illumination and detection geometry. Figures 6(b) and 6(c) show the images at $λ_{1}$ and $λ_{2}$ . The detector focused to a depth between the two pairs of blood vessels with a numerical aperture of 0.4 at a focal length of 2.5 cm. Electrical characteristics of the transducer and tissue attenuation were the same as in the example above. Table 1 lists the values of oxygenation and actual absorption coefficient. One approach to derive oxygenation levels would use ratio values of the pressure amplitudes at two wavelengths. Under the assumption of similar fluence, these values would be equal to the ratio of absorption coefficients. To assess the validity of this simplified approach in our simulation, we calculated the absorption coefficient ratios as well as the ratios of the pressure values after normalization with the incident radiant exposure. Signal amplitudes were obtained from areas of maximum amplitude in the B-scan images. At the chosen wavelengths, the ratio values agree quite well for the vessel with the higher oxygenation level. For the lower saturation level, however, the pressure amplitude ratio is lower than the corresponding absorption coefficient ratio. This is in agreement with the common observation that a simple ratio analysis that does not take into account the variation of local fluence with wavelength leads to erroneous oxygenation levels.³⁸

Fig. 6

Spectrally resolved PAM of blood vessels with different oxygen saturation ${sO}_{2}$ . (a) Outline of energy density distribution at 750 nm and sensor in transmission mode. (b) and (c) B-scan images at 797 and 750 nm, respectively.

Table 1

Spectral imaging of blood vessels with different oxygenation saturation levels sO2. Absorption coefficients μa,750 and μa,797 are given for the two wavelengths, 750 and 797 nm used in the simulations. The ratio of absorption coefficients is compared with the ratio of pressure amplitudes taken from the simulation in Fig. 6. z denotes the depth below the phantom surface, as shown in Fig. 6.

	sO2=95%	sO2=75%
$μ_{a, 750}$	$3.04 {cm}^{- 1}$	$3.99 {cm}^{- 1}$
$μ_{a, 797}$	$4.28 {cm}^{- 1}$	$4.28 {cm}^{- 1}$
$μ_{a, 750} / μ_{a, 797}$	0.71	0.93
$z = 2.4 c m$ : $p_{750} / p_{797}$	0.70	0.88
$z = 2 c m$ : $p_{750} / p_{797}$	0.72	0.86

3.3.

Comparison with Experiment

Figure 7 compares the measured and simulated B-scans of the tube phantom. In both cases, the same reconstruction using the transposed impulse response matrix was applied to the pressure signals. The ring-shaped conical detector is responsible for the “ $X$ ”-shaped appearance of the two absorbers in the images. At the focus line, the intersecting signals form the actual image of the object with a lateral resolution that only slowly changes with depth. This capability of generating a large depth of field is the benefit of this kind of axicon sensor. The drawback is the “ $X$ ”-shaped artifact, which can be reduced by some kind of deconvolution²⁵ or by using an array of concentric rings.²³^,³⁵ The width of the $X$ ’s, the depth dependence of amplitudes, and the overall appearance of the images agree between simulation and experiment.

Fig. 7

Comparison between simulated and experimental B-scan images of the phantom shown in Fig. 3. (a) Simulation and (b) experiment.

4. Discussion

The presented simulation method for photoacoustic imaging of an arbitrary distribution of deposited energy is particularly suited for situations, where the received acoustic waves have to be integrated over the surface of a finite size detector. In the specific examples of focused detector scanning, we could show that the simulation of a 2-D photoacoustic image from a Monte–Carlo generated absorbed energy distribution takes only about one minute. The spatiotemporal impulse response matrix and the energy density distribution have to be precomputed, requiring about one minute and some tens of minutes, respectively. The precomputed $H$ and $W (r)$ can be combined almost arbitrarily, meaning that one energy density distribution can be scanned with several sensors having different response matrices, or the same sensor can be tested with different $W (r)$ distributions. The only prerequisite is the use of appropriate discretization for the used acoustic frequency range. The system response used here is represented by a single matrix that models the complete process of wave propagation, spatial averaging over the sensor area and electrical characteristics of the transducer. Even an extension of this matrix taking into account a more complex behavior of the sensor, such as frequency-dependent directivity effects of a piezoelectric material would be possible. As the system response matrix is precomputed, such a modification would not affect the complexity of the actual simulation of pressure signals.

Modeling of photoacoustic signals from arbitrary sources is also possible with the pseudospectral method implemented in the k-wave toolbox for MATLAB. It is, therefore, worthwhile to compare the prerequisites and performance of the SIR method described in this study with $k$ -wave. First of all, there is a fundamental difference in how the signals are calculated: $k$ -wave computes the pressure field within the whole computational grid at defined time steps. In the time-domain integration underlying the SIR method, the signal at a defined receiving point is obtained for all times. To implement a finite size detector with $k$ -wave, e.g., a concave surface for focusing, the pressure is recorded for defined points lying on the sensor surface while the wave passes. In the SIR method, the single point signal is convolved with the spatially varying impulse response of the sensor to achieve signals for the focusing detector. The great advantage of the $k$ -wave method is that not only the source strength but also all other relevant material parameters can be defined point by point on the computational grid, such as the sound speed or the acoustic attenuation coefficient. The SIR method, on the other hand, assumes constant (average) values of these parameters. Concerning computation time, it has already been mentioned that a single $A$ -scan with the SIR method takes about the same time as a single time step with the $k$ -space method. Calculation of an A-scan with a focused detector, therefore, takes much longer with the latter, because a large number of time steps have to be computed for recording the acoustic wave passing through the points that form the detector surface. These points have to be within the computational grid, together with the initial pressure distribution from the MC simulation. This is not necessary when using the SIR method, where the detector response can be calculated also for remote source positions. Here, both the sensor surface and the source distributions can be simulated on a very fine grid.

We presented an approximation for a scanning transducer, where the absorbing structures within the simulated tissue phantom were shifted relative to the transducer and the energy density was obtained by multiplying the shifted absorption coefficient distribution with the static fluence distribution of the “empty” phantom. This phantom had appropriate values of background optical properties and contained the epidermis as a structure that did not change during scanning. This lead to a considerable acceleration of the computation, where the MC simulation was the main bottleneck. In our code, the MC simulation for $10^{6}$ photons required about half an hour. This is almost inhibitive for calculation of B-scan images, where on the order of one hundred such scans are necessary. However, faster implementations of MC codes are known, where it becomes feasible to calculate a new light distribution for each scanning position. Another alternative to multiple MC simulations could be methods that approximate the influence of absorbing structures on the fluence distribution, such as perturbation MC models.³⁹ Also for the diffusion approximation, perturbation methods are known, which are able to model the influence of inserted absorbing and scattering structures into a background light distribution.⁴⁰

5. Conclusions

The presented simulation method for PAM combines MC simulations of light transport with a time-domain method for acoustic propagation and detection, which is based on the SIR of a focusing sensor. The proposed method can be useful for accurate modeling of PAM, helping in the design of transducer and illumination systems. Furthermore, it provides an accurate forward model for any kind of model-based reconstruction in PAM.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Acknowledgments

This work has been supported by the Austrian Research Promotion Agency (FFG), Project No. 846559 and by the Austrian Science Fund (FWF), Project No. P 28032.

References

1.

K. Maslov, G. Stoica and L. V. Wang, “In vivo dark-field reflection-mode photoacoustic microscopy,” Opt. Lett., 30 (6), 625 –627 (2005). https://doi.org/10.1364/OL.30.000625 OPLEDP 0146-9592 Google Scholar

2.

M. Omar, J. Gateau and V. Ntziachristos, “Raster-scan optoacoustic mesoscopy in the 25–125 MHz range,” Opt. Lett., 38 (14), 2472 –2474 (2013). https://doi.org/10.1364/OL.38.002472 OPLEDP 0146-9592 Google Scholar

3.

J. Bauer-Marschallinger et al., “Fiber-optic annular detector array for large depth of field photoacoustic macroscopy,” Photoacoustics, 5 1 –9 (2017). https://doi.org/10.1016/j.pacs.2017.01.001 Google Scholar

4.

K. Maslov et al., “Optical-resolution photoacoustic microscopy for in vivo imaging of single capillaries,” Opt. Lett., 33 (9), 929 –931 (2008). https://doi.org/10.1364/OL.33.000929 OPLEDP 0146-9592 Google Scholar

5.

L. H. Wang, S. L. Jacques and L. Q. Zheng, “MCML—Monte–Carlo modeling of light transport in multilayered tissues,” Comput. Methods Programs Biomed., 47 (2), 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F CMPBEK 0169-2607 Google Scholar

6.

K. Sivasubramanian et al., “Optimizing light delivery through fiber bundle in photoacoustic imaging with clinical ultrasound system. Monte Carlo simulation and experimental validation,” J. Biomed. Opt., 22 (4), 041008 (2017). https://doi.org/10.1117/1.JBO.22.4.041008 JBOPFO 1083-3668 Google Scholar

7.

Z. Xie, L. V. Wang and H. F. Zhang, “Optical fluence distribution study in tissue in dark-field confocal photoacoustic microscopy using a modified Monte Carlo convolution method,” Appl. Opt., 48 (17), 3204 (2009). https://doi.org/10.1364/AO.48.003204 APOPAI 0003-6935 Google Scholar

8.

S. L. Jacques, “Coupling 3D Monte Carlo light transport in optically heterogeneous tissues to photoacoustic signal generation,” Photoacoustics, 2 (4), 137 –142 (2014). https://doi.org/10.1016/j.pacs.2014.09.001 Google Scholar

9.

S. A. Ermilov et al., “Development of laser optoacoustic and ultrasonic imaging system for breast cancer utilizing handheld array probes,” Proc. SPIE, 7717 717703 (2009). https://doi.org/10.1117/12.812192 PSISDG 0277-786X Google Scholar

10.

B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). https://doi.org/10.1117/1.3360308 JBOPFO 1083-3668 Google Scholar

11.

B. A. Kaplan et al., “Monte-Carlo-based inversion scheme for 3D quantitative photoacoustic tomography,” Proc. SPIE, 1006 100645J (2017). https://doi.org/10.1117/12.2251945 PSISDG 0277-786X Google Scholar

12.

M. Heijblom et al., “Appearance of breast cysts in planar geometry photoacoustic mammography using 1064-nm excitation,” J. Biomed. Opt., 18 (12), 126009 (2013). https://doi.org/10.1117/1.JBO.18.12.126009 JBOPFO 1083-3668 Google Scholar

13.

P. K. Upputuri et al., “Super-resolution photoacoustic microscopy using photonic nanojets. A simulation study,” J. Biomed. Opt., 19 (11), 116003 (2014). https://doi.org/10.1117/1.JBO.19.11.116003 JBOPFO 1083-3668 Google Scholar

14.

W. Liu et al., “Correcting the limited view in optical-resolution photoacoustic microscopy,” J. Biophotonics, 11 (2), e201700174 (2018). Google Scholar

15.

K. Wang et al., “An imaging model incorporating ultrasonic transducer properties for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 30 (2), 203 –214 (2011). https://doi.org/10.1109/TMI.2010.2072514 ITMID4 0278-0062 Google Scholar

16.

C. M. Á. Araque et al., “Model-based optoacoustic image reconstruction of large three-dimensional tomographic datasets acquired with an array of directional detectors,” IEEE Trans. Med. Imaging, 33 (2), 433 –443 (2014). https://doi.org/10.1109/TMI.2013.2286546 ITMID4 0278-0062 Google Scholar

17.

C. M. Á. Araque et al., “Model-based optoacoustic imaging using focused detector scanning,” Opt. Lett., 37 (19), 4080 –4082 (2012). https://doi.org/10.1364/OL.37.004080 OPLEDP 0146-9592 Google Scholar

18.

S. Karmakar, M. Roy and R. K. Saha, “Photoacoustic imaging of nanoparticle- containing cells using single-element focused transducer. A simulation study,” IEEE Trans. Ultrason. Ferroelect. Freq. Control, 62 (3), 463 –474 (2015). https://doi.org/10.1109/TUFFC.2014.006786 ITUCER 0885-3010 Google Scholar

19.

J. A. Jensen and N. B. Svendsen, “Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers,” IEEE Trans. Ultrason., Ferroelect., Freq. Control, 39 (2), 262 –267 (1992). https://doi.org/10.1109/58.139123 Google Scholar

20.

J. A. Jensen, Ultrasound Imaging and Its Modeling, 135 –166 Springer, Berlin (2002). Google Scholar

21.

S. L. Jacques, T. Li and S. A. Prahl, “3D Monte Carlo simulation of heterogeneous tissues,” http://omlc.org/software/mc/mcxyz/index.html Google Scholar

22.

S. L. Jacques, “Optical properties of biological tissues. A review,” Phys. Med. Biol., 58 (11), R37 –R61 (2013). https://doi.org/10.1088/0031-9155/58/11/R37 PHMBA7 0031-9155 Google Scholar

23.

K. Passler et al., “Scanning acoustic-photoacoustic microscopy using axicon transducers,” Biomed. Opt. Express, 1 (1), 318 –323 (2010). https://doi.org/10.1364/BOE.1.000318 BOEICL 2156-7085 Google Scholar

24.

R. G. M. Kolkman et al., “Photoacoustic imaging with a double-ring sensor featuring a narrow aperture,” J. Biomed. Opt., 9 1327 –1335 (2004). https://doi.org/10.1117/1.1805556 JBOPFO 1083-3668 Google Scholar

25.

S. Gratt et al., “Photoacoustic imaging using a conical axicon detector,” Proc. SPIE, 7371 73710W (2009). https://doi.org/10.1117/12.831739 PSISDG 0277-786X Google Scholar

26.

A. Penttinen and M. Luukkala, “The impulse response and pressure nearfield of a curved ultrasonic radiator,” J. Phys. D: Appl. Phys., 9 (10), 1547 –1557 (1976). https://doi.org/10.1088/0022-3727/9/10/020 JPAPBE 0022-3727 Google Scholar

27.

G. Paltauf and R. Nuster, “Artifact removal in photoacoustic section imaging by combining an integrating cylindrical detector with model-based reconstruction,” J. Biomed. Opt., 19 (2), 026014 (2014). https://doi.org/10.1117/1.JBO.19.2.026014 JBOPFO 1083-3668 Google Scholar

28.

K. Gurumurthy, “A dispersive model for the propagation of ultrasound in soft tissue,” Ultrasonic Imaging, 4 (4), 355 –377 (1982). https://doi.org/10.1016/0161-7346(82)90019-0 ULIMD4 0161-7346 Google Scholar

29.

G. F. Margrave, “Theory of nonstationary linear filtering in the Fourier domain with application to time‐variant filtering,” Geophysics, 63 (1), 244 –259 (1998). https://doi.org/10.1190/1.1444318 GPYSA7 0016-8033 Google Scholar

30.

B. E. Treeby, “Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering,” J. Biomed. Opt., 18 (3), 036008 (2013). https://doi.org/10.1117/1.JBO.18.3.036008 JBOPFO 1083-3668 Google Scholar

31.

R. Ma et al., “Fast scanning coaxial optoacoustic microscopy,” Biomed. Opt. Express, 3 (7), 1724 –1731 (2012). https://doi.org/10.1364/BOE.3.001724 BOEICL 2156-7085 Google Scholar

32.

E. Alerstam, T. Svensson and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt., 13 (6), 060504 (2008). https://doi.org/10.1117/1.3041496 JBOPFO 1083-3668 Google Scholar

33.

Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express, 17 (22), 20178 –20190 (2009). https://doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 Google Scholar

34.

L. Wang and S. L. Jacques, “Use of a laser beam with an oblique angle of incidence to measure the reduced scattering coefficient of a turbid medium,” Appl. Opt., 34 (13), 2362 –2366 (1995). https://doi.org/10.1364/AO.34.002362 APOPAI 0003-6935 Google Scholar

35.

G. Paltauf et al., “Ring detector arrays for large depth of field scanning photoacoustic macroscopy,” Proc. SPIE, 10494 104943Q (2018). https://doi.org/10.1117/12.2292051 PSISDG 0277-786X Google Scholar

36.

B. T. Cox and P. C. Beard, “The frequency-dependent directivity of a planar Fabry–Perot polymer film ultrasound sensor,” IEEE Trans. Ultrason. Ferroelect. Freq. Control, 54 (2), 394 –404 (2007). https://doi.org/10.1109/TUFFC.2007.253 ITUCER 0885-3010 Google Scholar

37.

V. A. Shutilov, Fundamental Physics of Ultrasound, Gordon and Breach, New York (1988). Google Scholar

38.

B. Cox et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt., 17 (6), 061202 (2012). https://doi.org/10.1117/1.JBO.17.6.061202 JBOPFO 1083-3668 Google Scholar

39.

A. Sassaroli, “Fast perturbation Monte Carlo method for photon migration in heterogeneous turbid media,” Opt. Lett., 36 (11), 2095 –2097 (2011). https://doi.org/10.1364/OL.36.002095 OPLEDP 0146-9592 Google Scholar

40.

M. R. Ostermeyer and S. L. Jacques, “Perturbation theory for diffuse light transport in complex biological tissues,” J. Opt. Soc. Am. A, 14 (1), 255 (1997). https://doi.org/10.1364/JOSAA.14.000255 JOAOD6 0740-3232 Google Scholar

Biography

Guenther Paltauf received his PhD in physics from the University of Graz, Austria, in 1993. After postdoctoral studies in Graz, the University of Bern, Switzerland and the Oregon Medical Laser Center in Portland, Oregon, he is now an associate professor at the Department of Physics at the University of Graz. His research interests are photoacoustic imaging, photomechanical laser-tissue interactions, and laser-ultrasonic characterization of materials.

Paul R. Torke completed his master’s degree in biomedical engineering from the Technical University Ilmenau, Germany, in 2012. From 2012 to 2015, he worked as a clinical specialist in the field of cardiological interventions. Currently, he is working on his PhD in experimental physics at the Karl-Franzens-University Graz, Austria.

Robert Nuster received his PhD in experimental physics from Karl-Franzens-University Graz, Austria, in 2007, with a thesis on development and application of optical sensors for laser induced ultrasound detection. From 2008 to 2011, he was a postdoctoral research fellow at the Karl-Franzens-University Graz. Since 2011, he has been a senior postdoctoral research fellow. His current research interests include photoacoustic imaging, characterization of materials by laser ultrasound, and ultrasound sensor development.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Guenther Paltauf, Paul R. Torke, and Robert Nuster "Modeling photoacoustic imaging with a scanning focused detector using Monte Carlo simulation of energy deposition," Journal of Biomedical Optics 23(12), 121607 (24 September 2018). https://doi.org/10.1117/1.JBO.23.12.121607

Received: 19 June 2018; Accepted: 15 August 2018; Published: 24 September 2018

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 16 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Sensors

Acoustics

Monte Carlo methods

Signal detection

Photoacoustic imaging

Photoacoustic spectroscopy

Signal attenuation

1.

Introduction

2.

Methods

2.1.

Monte Carlo Simulation of Light Transport

Eq. (1)

2.2.

Photoacoustic Sound Generation

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Fig. 1

2.3.

Scanning the Detector

Fig. 2

2.4.

Experiment

Fig. 3

3.

Results

3.1.

Illumination Modes

Fig. 4

Fig. 5

3.2.

Spectrally Resolved Imaging

Fig. 6

Table 1

3.3.

Comparison with Experiment

Fig. 7

4.

Discussion

5.

Conclusions

Disclosures

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years