Modeling toolchain for realistic simulation of photoacoustic data acquisition

Jan-Willem Muller; Mustafa Umit Arabul; Hans-Martin Schwab; Marcel C. M. Rutten; Marc R. H. M. van Sambeek; Min Wu; Richard G. P. Lopata

doi:10.1117/1.JBO.27.9.096005

14 September 2022 Modeling toolchain for realistic simulation of photoacoustic data acquisition

Jan-Willem Muller, Mustafa Umit Arabul, Hans-Martin Schwab, Marcel C. M. Rutten, Marc R. H. M. van Sambeek, Min Wu, Richard G. P. Lopata

Author Affiliations +

Journal of Biomedical Optics, Vol. 27, Issue 9, 096005 (September 2022). https://doi.org/10.1117/1.JBO.27.9.096005

Abstract

Significance: Physics-based simulations of photoacoustic (PA) signals are used to validate new methods, to characterize PA setups and to generate training datasets for machine learning. However, a thoroughly validated PA simulation toolchain that can simulate realistic images is still lacking.

Aim: A quantitative toolchain was developed to model PA image acquisition in complex tissues, by simulating both the optical fluence and the acoustic wave propagation.

Approach: Sampling techniques were developed to decrease artifacts in acoustic simulations. The performance of the simulations was analyzed by measuring the point spread function (PSF) and using a rotatable three-channel phantom, filled with cholesterol, a human carotid plaque sample, and porcine blood. Ex vivo human plaque samples were simulated to validate the methods in more complex tissues.

Results: The sampling techniques could enhance the quality of the simulated PA images effectively. The resolution and intensity of the PSF in the turbid medium matched the experimental data well. Overall, the appearance, signal-to-noise ratio and speckle of the images could be simulated accurately.

Conclusions: A PA toolchain was developed and validated, and the results indicate a great potential of PA simulations in more complex and heterogeneous media.

1. Introduction

Photoacoustic (PA) imaging is a relatively new biomedical imaging modality that employs the wide variety of optical properties found in tissue to generate contrast with ultrasound (US) resolution. Increasingly sophisticated PA setups and signal processing methods are being developed to enhance the capabilities of PA imaging. Phantom experiments are often performed to validate new methods, to characterize PA setups, and to generate training datasets for machine learning. However, the scarcity of high-quality phantoms and the absence of a ground truth in in vitro experiments impedes research. Numerical modeling approaches can provide an in silico research environment, in which all relevant physical parameters are known and controllable. Therefore, an accurate and validated PA imaging toolchain that can generate realistic PA images may accelerate the general development of PA imaging. However, generating realistic PA image data with simulations remains a great challenge. In this research, we propose a thoroughly validated, quantitative, and flexible PA imaging modeling toolchain and demonstrate its capability to simulate the acquisition of PA signals of a physical PA system in both phantoms and complex tissues.

The fundamental optical and acoustic processes underlying PA imaging are generally well-understood, allowing for physics-based modeling of PA signal generation and acquisition. As this kind of modeling is strictly governed by physical relations, it is inherently transparent and enables full control over the simulated phenomena. Physics-based PA modeling techniques have been used for a wide variety of purposes, e.g., to investigate model-based image reconstruction, to optimize PA setups, and to analyze acoustic effects and artifacts, such as speed of sound aberrations and clutter.¹^–⁵ Another emerging research area in which the simulation of data plays a major role is the development of machine learning to improve or interpret measured PA data. Machine learning has the potential to increase the quality of PA images by supressing clutter and noise, to overcome limited-view and limited-bandwidth artefacts, and to provide quantitative estimates of the tissue’s optical and acoustic properties.⁶ However, machine learning is a data-driven approach, and a major bottleneck in its application is the lack of reliable experimental training data and the ground truth. The use of accurate numerical models is therefore essential to generate reliable training datasets in machine learning.⁷

In biomedical applications, the duration of the optical pulses used is typically short enough to meet the thermal and stress confinement criteria needed for adequate optical energy delivery and image resolution.⁸ When these criteria are met, the modeling of PA signals is generally subdivided into an optical part and an acoustic part. First, the locally resolved pressure rise, caused by the absorption of laser light and thermoelastic expansion of the tissue, is determined. To calculate the initial pressure, the optical fluence, which is physically governed by the radiative transfer equation (RTE), needs to be determined. Solving the RTE in light–tissue interaction is done either with statistical Monte-Carlo (MC) methods or using a diffusion approximation of the RTE, which can be solved with continuum methods (e.g., finite-difference and finite element methods).⁹^–¹¹ Next, the propagation of the pressure waves must be simulated by solving the acoustic wave equations. Many different mathematical techniques and software packages to model wave propagation of medical US have been described in the literature, often relying on spectral methods or the finite element method.¹² For efficient and accurate PA wave propagation modeling and reconstruction in heterogeneous media, the open source toolbox k-Wave is frequently used in the literature.⁷^,¹³^,¹⁴

To be used as a reliable tool to generate realistic PA data, a PA imaging toolchain must be flexible enough to adapt to various setups and media. Typical PA setups can consist of a range of laser sources, e.g., fiber lasers,¹⁵ crystal lasers,¹⁶ and diode lasers,¹⁷^,¹⁸ and different types of the acoustic sensors, including single element transducers,¹⁹ linear array transducers,¹⁸ and custom-designed transducers.²⁰^,²¹ Furthermore, the tissues on which PA imaging is performed, are often heterogeneous and induce both optical and acoustic scattering, reflections, refraction, and/or attenuation.

The coupling of the optical and acoustic domains, as well as a correct description of both the PA setup and the medium to be imaged are nontrivial tasks. A thorough analysis of the quality of the simulated PA data is needed to ensure realistic results. These analyses, however, are often qualitative or rely on idealized setups and/or simplified media, such as wire phantoms.⁴^,⁷^,¹⁴^,²²^–²⁵ A thorough validation of the entire modeling chain with physical PA set-ups in complex media is still missing. This lack of validation makes the accuracy and the realism of the simulations an uncertainty and thus hinders the translation toward in vivo applications.

Therefore, in this study, a toolchain is developed and validated by comparing simulated PA images with those acquired by an experimental PA imaging setup. To further improve the quality of the simulated images, this study also investigates the performance of a sampling strategy for the definition of the initial pressure field prescription in k-Wave. The toolchain is designed to be quantitative, by simulating channel data in an absolute unit (Pascal) with a physical model. Furthermore, it was made to be as flexible as possible with respect to the medium geometry and PA system specifications. Both the optical and acoustic parts of the toolchain are validated, and its capability to simulate realistic PA images is demonstrated using multiple in vitro phantoms and ex vivo human tissue samples of carotid plaques.²⁶

2. Methods

2.1.

Optical Fluence Modeling

A grid-based MC software was developed to simulate the optical fluence distribution in complex media. The tool is implemented in C++, callable from MATLAB (R2021a, The Mathworks, Natick, Massachusetts) and compiled with multithreading support using OpenMP (GCC 8.1, Free Software Foundation, Boston, Massachusetts). The implementation is based on the microscopic Beer–Lambert law method, which shows an overall superior convergence rate compared with other methods.²⁷^,²⁸ In this method, photon packets with energy $E_{p}$ are tracked that deposit energy in the medium in a continuous way:

Eq. (1)

E_{d} = E_{p} \cdot (1 - e^{- μ_{a} \cdot L_{a}}),

where

E_{d}

is the energy (J) deposited to a voxel in the medium,

L_{a}

is the distance (m) traveled by the photon packet in that voxel, and

μ_{a}

is the optical absorption coefficient (

m^{- 1}

). Energy conservation is maintained by subtracting

E_{d}

from

E_{p}

after each deposition. The distance traveled between two scattering events,

L_{s}

(m), is determined as follows:

Eq. (2)

L_{s} = - \frac{\log (η)}{μ_{s}},

where

μ_{s}

denotes the scattering coefficient (

m^{- 1}

) and

η

is a random number generated from a uniform probability distribution between 0 and 1. The direction of the photon packet is changed in a scattering event according to the Henyey–Greenstein phase function:

Eq. (3)

p (θ) = \frac{(1 - g^{2}) \cdot \sin (θ)}{2 \cdot (1 + g^{2} - 2 \cdot g \cdot \cos (θ))^{\frac{3}{2}}},

where

g

is the scattering anisotropy (–),

θ

is the angular deflection with respect to the photon packet’s previous direction (rad), and

p (θ)

is its corresponding probability function.

The optical properties for each grid point in the domain ( $μ_{a}$ , $μ_{s}$ , $g$ ) can be set according to the local simulated tissue compositions, allowing for a simulation of the photon distribution in complex heterogeneous media.

2.2.

Acoustic Simulation and Sensor Definition

The propagation of the initial pressure field in the tissue over time is simulated according to the linear wave equation using the k-Wave toolbox. k-Wave is an open source toolbox that uses a pseudospectral, time-domain, finite difference method to solve the linear wave equation for compressional waves.¹³ The largest supported frequency in the simulation is limited by the Nyquist criterion, which states that a sampling of at least two points per wavelength (PPWL) is needed. The largest allowed spacing in the simulation is therefore related to the maximum frequency that the US probe can detect. All waves with a higher frequency will not be sensed by the probe and can thus be safely neglected in the simulation.

To define a sensor in k-Wave, each transducer element must be discretized to pixel positions and associated weights. The position of an element is usually mapped to the grid using nearest-neighbor interpolation.¹³ This approach, however, alters the size and position of the elements and the kerf between them. Furthermore, it causes staircasing artifacts when the elements are not aligned with the grid, negatively affecting the simulated channel data. These artifacts can be reduced by band-limiting the shape of the elements first, using a numerical convolution with the band-limited interpolant (BLI), which gives the correct quadrature weights for discretization.²⁹ In this study, an analytical description of the band-limited rectangular element shape was used. The shape was prescribed directly in the spatial Fourier domain using a sinc function, which gave the appropriate quadrature weights after applying the inverse Fourier transform.

2.3.

Initial Pressure Definition

The Nyquist criterion needs to be considered when sampling $p_{0}$ for use in k-Wave. If the maximum spatial frequency of the initial pressure is higher than the Nyquist frequency, aliasing occurs when sampling this distribution. For initial pressure distributions of structures not aligned with the grid, this effect may appear as staircasing artifacts (pixelated/discontinuous edges) in the reconstructed image.

To avoid aliasing errors, an analytical description of the initial pressure in the spatial frequency domain may be used. The initial pressure field can subsequently be sampled in the frequency domain, which prevents the occurrence of aliasing. This approach, however, is limited to spatial pressure distributions for which analytical descriptions exist.

A more generic and straightforward method would be to use a finer spatial sampling of $p_{0}$ and adjust the acoustic simulation grid accordingly. This method will be referred to as the direct simulation method. Although a direct simulation will converge to the correct solution for decreasing $Δ x_{sim}$ , it is not viable for most purposes, as it increases computation time significantly, especially for two-dimensional (2D) and three-dimensional (3D) simulations.

Wise et al.²⁹ introduced a technique in which the band-limited pressure field is obtained by a numerical integration of point sources, using the BLI. Although this method is theoretically valid for all geometrical shapes, finding the optimal spatial distribution of point sources for arbitrary shapes is a nontrivial task. Furthermore, the BLI has an infinite support, which makes this method computationally demanding. To reduce computation time, truncation of the BLI was suggested at the expense of accuracy.

In this study, we propose a computationally inexpensive and accurate preprocessing step to reduce aliasing artifacts, see Fig. 1. The initial pressure (in Pa) is calculated according to

Eq. (4)

p_{0} = Γ \cdot μ_{a} \cdot Φ,

where Γ denotes the Grüneisen parameter [–],

μ_{a}

is the optical absorption coefficient (

m^{- 1}

) and Φ is the optical fluence (

J / m^{2}

) as determined by the optical MC simulation. Unlike the relatively large spatial spacing used for the acoustic simulation,

Δ x_{sim}

, a smaller spacing spatial spacing

Δ x_{pre}

(

Δ x_{pre} < Δ x_{sim}

) is used to sample

p_{0}

. This smaller

Δ x_{pre}

allows for a better sampling of

p_{0}

, as less content in the spatial frequency domain is aliased.

Fig. 1

A schematic example of the proposed $p_{0}$ sampling strategy of a rectangular absorber. After determining the optical fluence, $p_{0}$ is calculated using a fine grid spacing. A filtering and resampling step is performed subsequently to sample $p_{0 sim}$ on the coarse grid used in the acoustic simulation.

To create speckle in the PA images, a Gaussian distribution was used to model deviations in the absorption coefficient and Grüneisen parameter. In this model, random fluctuations were added to the initial pressure:

Eq. (5)

p_{0_{speckle}} = p_{0} + N (σ_{speckle}) \cdot Φ,

Eq. (6)

σ_{speckle} = \sqrt{I_{speckle}} \cdot (Δ x_{pre})^{- \frac{D}{2}},

where

N (σ)

denotes a zero mean normal distribution with standard deviation

σ

,

I_{speckle}

is a spatially distributed quantity proportional to the local acoustic speckle intensity, and

D

is the number of simulated dimensions (

D = 1, 2, 3

, corresponding to 1D, 2D, and 3D simulations) in the acoustic simulation. The correction factor

{(Δ x_{pre})}^{- (D / 2)}

ensures the simulated speckle intensity is independent of the used sampling spacing, see Appendix A.

Next, $p_{0_{speckle}}$ is filtered with a multidimensional ideal low-pass filter $H$ to remove frequencies higher than the maximal frequency supported in the simulation:

Eq. (7)

rect (ξ) = {\begin{cases} 1 & for | ξ | \leq \frac{1}{2} \\ 0 & for | ξ | > \frac{1}{2} \end{cases},

Eq. (8)

H (\vec{k}) = \prod_{i = 1}^{D} rect (k_{i} \cdot Δ x_{sim}),

Eq. (9)

p_{0_{speckle}}^{*} = F^{- 1} {F (p_{0_{speckle}}) \cdot H (\vec{k})},

where

F

denotes the multidimensional spatial Fourier transform and

\vec{k}

spatial frequency. Finally, the filtered distribution

p_{0_{speckle}}^{*}

is downsampled to obtain

p_{0 sim}

, which was prescribed as input in k-Wave. A schematic overview of all processing steps in the toolchain, is shown in Fig. 2.

Fig. 2

A schematic overview of all processing steps. The gray boxes indicate enhanced processing methods proposed in this study.

2.4.

Staircasing Error Reduction Analysis

To estimate the performance of the method proposed, both the required computation time and the error in the delay-and-sum reconstructed PA images of disc-shaped $p_{0}$ distributions, with a radius of 1 and 5 mm, were determined. The simulations were performed on a PC with an i7-10700 processor (Intel, Santa Clara, California) and an RTX 3070 GPU (Nvidia, Santa Clara, California), using the optimized CUDA k-Wave code. The $p_{0 sim}$ distributions were obtained using three sampling strategies: (1) the preprocessing method proposed with a fixed $Δ x_{sim}$ and a varying $Δ x_{pre}$ , (2) by means of the Wise et al. method²⁹ with a fixed $Δ x_{sim}$ and a varying number of integration points $N_{int}$ , and (3) using a direct simulation with varying $Δ x_{sim}$ . The kWaveArray class (alpha version 0.3) was used as implementation of the Wise et al. method, and the default truncation of the BLI at 5% was applied. In this implementation, the disc is sampled in a radial pattern, with the number of integration points increasing linearly as a function of the radial position to enhance uniformity. For comparison, $N_{int}$ was expressed as a function of the effective sampling distance $Δ x_{eff}$ using the following equation:

Eq. (10)

N_{int} = \frac{A_{obj}}{(Δ x_{eff})^{2}},

where

A_{obj}

denotes the area (

m^{2}

) of the simulated object.

The simulation results were compared with an analytical solution of a low-pass filtered disc phantom that exhibits no staircasing. This reference PA phantom was obtained by means of a simulation in which $p_{0 sim}$ was sampled in the spatial frequency domain:

Eq. (11)

p_{0_{sim}} = F^{- 1} (\frac{J_{1} (2 π \sqrt{ξ_{x}^{2} + ξ_{y}^{2}})}{\sqrt{ξ_{x}^{2} + ξ_{y}^{2}}}),

where

ξ_{x}

and

ξ_{y}

denote the spatial frequency (scaled with respect to the radius) in

x

and

y

directions, respectively, and

J_{1}

is the order-1 Bessel function of the first kind.

The error $ε$ , with respect to the reference simulation, was calculated using

Eq. (12)

ε = 10 \cdot \log_{10} \frac{\sum (p - p_{ref})^{2}}{\sum p_{ref}^{2}},

where

p

and

p_{ref}

denote the pressures of the reconstructed radio frequency signals, respectively.

2.5.

In Vitro and In Silico Experimental Setup

The developed toolchain was validated with a real experimental setup. The setup consisted of a tunable pulsed laser (Opotek, Radiant HE 355LD, Carlsbad, California) with a pulse duration of 6 ns and a fiber bundle with a custom designed circular output aperture of 4.9 mm in diameter (CeramOptec, Bonn, Germany). A Vantage 256 system (Verasonics, Kirkland, Washington) connected to a Verasonics L11-5v linear array transducer with 128 elements, a pitch of $300 μ m$ , and an elevational focus depth of 18 mm, was used as US detector. The probe’s center frequency was 7.6 MHz, and it had a bandwidth of 75%. The frequency response of the probe was implemented in the toolchain using a linear frequency filter, defined as a Hann window in the temporal Fourier domain. The center frequency and bandwidth of the Hann filter (at $- 6 dB$ ) in the simulations were chosen to match the probe used. The pixel spacing used in k-Wave, $Δ x_{sim}$ , was set to $30 μ m$ , unless stated otherwise. The fiber bundle and probe were aligned and mounted to a linear stage using a custom-designed probe holder, see Figs. 3(a) and 3(b).³⁰ In the optical simulations, the fiber bundle tip was approximated as a disc source with uniform intensity. The angular distribution of the emitted photons was uniform within a cone with a maximum angle with respect to the bundle’s surface of 5 deg, see Figs. 3(c) and 3(d). All optical simulations were performed in 3D with a grid spacing of $250 μ m$ . The medium properties were constant in the elevational direction. The acoustic simulations, however, were performed in 2D in the imaging plane of the linear probe. Acoustic attenuation was not considered, as the acoustic attenuation is relatively insignificant compared with the effective optical attenuation for all experiments.

Fig. 3

(a) A schematic drawing of the setup used. The fiber bundle and probe are attached to a linear stage using the probe holder. The mounted sample can be rotated by a motor. (b) A photograph of the probe holder. (c) A simulation of the normalized fluence ( $m^{- 2}$ ) in water. The probe’s imaging plane is indicated by the dashed line. (d) The optical fluence in the probe’s imaging plane, corresponding to the green dashed line in (c).

2.6.

Verification of Simulated PSF in Turbid Media

The ability of the toolchain to accurately simulate the point spread function (PSF) in turbid media was verified experimentally. A black thread ( $280 μ m$ diameter) served as an absorber and was positioned in a tank filled with deionized water. The angle between the normal direction of the probe and the laser source (800 nm) was 40 deg. The position of the probe and the laser source could be adjusted simultaneously using a linear stage connected to the probe holder, to vary the distance to the black thread in steps of 0.5 mm, see Fig. 3. The reconstructed pressure at the thread’s position in the PA images was determined for each position and were compared between the in vitro and simulated results.

Intralipid (Intralipid 20%, Fresenius Kabi Nederland, Huis Ter Heide, The Netherlands) was added to the water to increase the scattering coefficient. The 20% intralipid was diluted to solutions with a range of five concentrations (0.01%, 0.02%, 0.05%, 0.1%, and 2%). This range of concentrations corresponds to reduced scattering coefficients between 0.1 and $2 {cm}^{- 1}$ .³¹ As the absorption of intralipid at 800 nm wavelength is negligible, only absorption due to the water ( $μ_{a} = 2.3 m^{- 1}$ ) was taken into account in the simulations.³¹^,³² The scattering anisotropy factor was set to 0.55, in agreement with estimations of the anisotropy factor of intralipid using Mie theory³¹.

2.7.

Three-Channel Phantom

To assess the performance of the acoustic simulation in reproducing the visual appearance of realistic tissues, a controlled in vitro study was conducted using a polyvinyl alcohol (PVA) phantom with three channels representing plaque in a carotid artery, see Figs. 4(a)–4(d).³³ The three channels were filled with: (1) cholesteryl linoleate (Sigma-Aldrich, St. Louis, Missouri), (2) pieces of a human plaque, and (3) porcine blood collected from the slaughterhouse with 19% w/v sodium citrate tribasic dihydrate (Sigma-Aldrich, St. Louis, Missouri) as anticoagulant. The human plaque samples in this study were obtained during a carotid endarterectomy as part of an in vivo PA imaging study previously reported by our group.²⁶ The phantom was fixed in a rotatable set-up that allowed for tomographic compounding to increase PAI quality and field of view.³³ An optical wavelength of 710 nm was chosen to illuminate the phantom, as this resulted in relatively high-quality images with uniform absorption within each channel. The phantom was rotated over 180 degrees in steps of 18 degrees, to apply an even illumination. A uniform fluence distribution was used in the acoustic simulations. Gaussian noise was added to the simulated channel data to match the noise level of the in vitro measurements and the Grüneisen parameter was spatially homogeneous and fixed for all simulations. An identical band-pass filter (3 to 13 MHz) was applied to both the simulated and the in vitro channel data to remove noise and artifacts beyond the bandwidth of the transducer. Afterward, the channel data were delay-and-sum beamformed to create an image for each rotation angle. The resulting images were spatially compounded, either coherently (averaging of radio frequency data) or incoherently (averaging after envelope detection).³⁴ To mimic registration errors in the in vitro set-up, the simulated PA images were translated randomly in both axial and lateral direction, using an empirically determined Gaussian distribution with a standard deviation of $200 μ m$ before spatial compounding.

Fig. 4

The properties of the three-channel phantom used in the toolchain validation. (a) The three channels are filled with cholesterol, a human plaque sample, and blood, respectively. (b) The initial pressure field, (c) the speed of sound map, and (d) the density map as prescribed in the acoustic simulation. (e) and (f) Two in vitro US images of human plaque samples. The corresponding segmentations used in the simulations are shown in (g) and (h).

2.8.

Ex Vivo Plaque Imaging

The developed toolchain was used to simulate ex vivo PA images of two human plaque samples to evaluate the performance in a more complex and heterogeneous tissue. The plaques were scanned using an optical wavelength of 800 nm. Manual segmentations of the plaque were made using US B-mode data and histology. Different tissue types were assigned to each region in the plaque, see Figs. 4(e)–4(h). The optical and acoustic properties in the simulation were assigned per region, see Table 1. Note that the tissue types and corresponding material properties for each region were tuned to match the ex vivo PA images and may not reflect the ground truth exactly. The. $I_{speckle}$ and the noise level were tuned empirically to be in agreement with the observed speckle intensity and noise in the ex vivo images and the Grüneisen parameter was spatially homogeneous and fixed for both simulations. An identical band-pass filter (3 to 13 MHz) was applied to both the simulated and the in vitro channel data.

Table 1

The optical (at 800 nm wavelength) and acoustic parameters used in each segmented region of the two plaques. The reference value for Ispeckle at 0 dB is 8.9·10−6.

Tissue type	Speed of sound (m/s)	Density (kg/m3)	Ispeckle (dB)	Optical absorption (1/cm)	Optical reduced scattering (1/cm)	Anisotropy factor g (∼)
Lumen/water (1)	1480	1000	$- \infty$	0.023³²	0³²	—
Fibrous tissue (2)	1540	1050	$- 10$	0.8³⁵	12.6³⁶	0.93³⁶
Necrotic core (3)	1540	1000	$- 10$	0.8³⁵	13.4³⁶	0.93³⁶
Hemorrhage (4)	1540	1000	$- 5$	1.6³⁵	12.6³⁶	0.93³⁶
Fibrous tissue (5)	1540	1050	0	0.8³⁵	12.6³⁶	0.93³⁶
Calcification (6)	1580	1100	0	1.1³⁵	16.4³⁶	0.93³⁶

3. Results

3.1.

Initial Pressure Preprocessing

A simulation of the 5 mm radius disc with the default spatial sampling width $Δ x_{sim}$ of $30 μ m$ (6.6 PPWL at the probe’s center frequency) leads to significant aliasing artefacts in the reconstructed image, see Fig. 5(a). Although the top and bottom signals are reconstructed correctly, staircasing artifacts arise in the regions between 1 and 5 o’clock, and between 7 and 11 o’clock. Furthermore, the side lobes and grating lobes of these staircased signals cause an erroneously elevated background signal in a large part of the image. By increasing the PPWL to 40, which results in a $p_{0}$ spatial sampling width of $4.9 μ m$ , the errors could be reduced significantly, at the cost of a strong increase in memory occupation and computation time. The computation time increased from $\sim 14 s$ to 36 min between the simulations performed with 6.6 PPWL and 40 PPWL using the direct method. The method by Wise et al.²⁹ increased the computation time from 22 s to 4.7 min and showed an erroneously elevated signal within the disc, but no staircasing artifacts at the disc’s boundary. The PA image simulated using the preprocessing method is almost identical to the one simulated with the direct method. Hence, the preprocessing step is effectively reducing aliasing errors. The additional computation time needed for the preprocessing step was 3 s, yielding a total simulation time of 17 s.

Fig. 5

The reconstructed PA images of the 5-mm radius disc are shown (60 dB dynamic range) for different sampling strategies of $p_{0}$ : (a) $p_{0}$ was sampled directly on the k-Wave grid with $30 μ m$ spacing (6.6 PPWL); (b) $p_{0}$ was directly sampled on a finer k-Wave simulation grid with 40 PPWL; (c) $p_{0}$ was sampled using the preprocessing method proposed with 40 PPWL; (d) $p_{0}$ was sampled using the Wise et al. method with 40 PPWL. The reference PA image is shown in (e), in which $p_{0}$ was sampled in the spatial Fourier domain.

The effectiveness of both the direct method and the preprocessing method was quantified by calculating the error of the PA images, $ε$ , with respect to the results from the reference simulation, see Fig. 6. The reduction of the error is similar between the direct simulation and the preprocessing method, for an equal number of PPWL. For these methods, $ε$ is independent of the size of the disc-shaped $p_{0}$ distribution and the center frequency of the probe and shows algebraic convergence with increasing PPWL. For the method proposed by Wise et al.,²⁹ the results show a similar trend, although the intensity error is less similar between different disc sizes and center frequencies, see Fig. 6(c). The simulations using the preprocessing method proposed were generally faster than those using the Wise et al. method, up to 47 times for the large disc at 7.6 MHz center frequency with 80 PPWL, see Fig. 6(d). An equation was fitted empirically to determine the mean error of the reconstructed image as function of the number of PPWL for the preprocessing method proposed:

Eq. (13)

ε \approx - 13.3 \cdot \ln (PPWL) + 9.09 .

Fig. 6

(a) Intensity error $ε$ is shown as a function of the pixel spacing of the acoustic simulation using the direct simulation method. (b) $ε$ is shown as function of $Δ x_{pre}$ in the method proposed,²⁹ with a constant acoustic simulation pixel spacing of $30 μ m$ . (c) ε is shown as a function of $Δ x_{eff}$ in the method by Wise et al., with a constant acoustic simulation pixel spacing of $30 μ m$ . The solid red lines show the fit as described by Eq. (13). (d) The computation time ratio needed for the combined $p_{0}$ sampling and k-Wave simulation between the method by Wise et al.²⁹ and the method proposed.

3.2.

Fluence Simulation in Turbid Medium

The PSF of the black thread was measured in vitro and compared with the simulated result, see Figs. 7(c) and 7(d). Overall, the appearances of the in vitro and simulated PSF are similar. The axial resolution of the in vitro PSF, $310 μ m$ ( $- 6 dB$ peak width) was slightly worse than the simulated resolution of $270 μ m$ . The in vitro and simulated positions of the grating lobes are almost identical, and the lobes’ amplitudes are comparable in general. The amplitude of the reconstructed pressure was measured for a range of intralipid concentrations in vitro, see Figs. 7(a) and 7(b). By increasing the intralipid concentration, the reduced scattering coefficient is increased as well, leading to a smaller peak amplitude and a shift of the peak location toward smaller depths. The change in the peak’s amplitude and depth as function of intralipid concentration in the simulated data corresponded well to the experimental results over a relatively large range of pressure amplitudes (30 dB).

Fig. 7

(a) The maximum reconstructed pressure of the black thread phantom is shown for water and two different concentrations of intralipid. (b) The maximum pressures and corresponding depths of the PSF are shown for water and all concentrations of intralipid. The dashed line shows the results for additionally simulated concentrations of intralipid. (c) The simulated PSF in water. (d) The in vitro PSF in water. A dynamic range of 60 dB was used for visualization.

3.3.

Three-Channel Phantom

The in vitro and simulated PA images of the three-channel phantom were obtained by spatial compounding of the acquisitions over all rotation angels, see Figs. 8(a)–8(d). The in vitro and simulated PA images appear to be similar in terms of overall contrast and resolution of the PVA phantom and its three channels for both compounding techniques. The signal-to-noise ratio (SNR) was determined for all channels and compounding techniques, see Figs. 8(e) and 8(f). Overall, the SNR of the coherent compounded simulated images in which no registration errors were added was more than 5 dB larger than the in vitro data. The increased SNR was observed in the incoherent simulated images too, although less pronounced. By including registration errors in the simulated data, the mean SNR of all channels was decreased by $4.4 \pm 0.39 dB$ for coherent compounding, and by $0.77 \pm 0.10 dB$ , compared with simulated data with perfect registration. After adding registration errors, the SNR of simulated data agreed well with the in vitro data for all channels and for both compounding techniques. The mean SNR decrease for coherent compounding was 3.6 dB larger than the decrease for incoherent compounding. This observation can be explained by the fact that coherent averaging is more susceptible to registration errors than incoherent compounding, as the latter does not take phase information into account.

Fig. 8

(a), (b) The simulated PA images for incoherent and coherent spatial compounding. (c), (d) The in vitro PA images for incoherent and coherent spatial compounding. All images are shown with a dynamic range of 30 dB. The SNR ( $n = 10$ for simulated data, error bars denote standard deviation) for each channel is shown in (e) for the coherently compounded images, and in (f) for the incoherently compounded images. The SNR of the simulated data was calculated both with an imperfect registration (I.R.) and with a perfect registration (P.R.).

3.4.

Ex Vivo Plaque Imaging

Two human carotid plaques were scanned ex vivo at a wavelength of 800 nm, see Figs. 9(a) and 9(b). The top parts of the plaques show relatively strong speckle signals. The intensity of these signals drops rapidly over only several millimetres deeper into the plaque, implying a strong fluence decrease with increasing depth. Furthermore, two PA interface signals originating from sharp transitions in the local absorption coefficient can be observed in the first plaque, see Fig. 9(a). These signals, indicated by the green arrows, appear only as lines due to the band-limitation of the probe.²⁶ Overall, the simulated PA images agree well with the ex vivo images in terms of overall appearance, resolution, and contrast, see Figs. 9(c) and 9(d). Furthermore, both the rapid fluence drop as function of depth, and the two PA interface signals could be simulated realistically. However, the top region of the plaque in Fig. 9(b) shows several high-intensity absorbers that were not present in the used segmentation. The Gaussian model for speckle could not simulate these high-intensity absorbers.

Fig. 9

(a), (b) Two ex vivo PA images (at 800 nm) of different human plaques (P1, P2). (c), (d) The corresponding simulations. Two PA signals originating from the interface between two absorbing layers can be observed in (a) and (c) and are indicated by the arrows. The dynamic range of all images is 35 dB.

4. Discussion

Many numerical methods to simulate the acoustic wave propagation and the optical absorption in tissue have been described and verified before, using both analytical solutions and in vitro experiments.¹¹^,¹²^,³⁷^,³⁸ However, a thorough validation of the entire modeling chain with physical PA setups in complex media is still missing. This lack of validation makes the accuracy and the realism of the simulations an uncertainty and thus hinders the translation toward in vivo applications. In this study, a modeling toolchain has been developed and thoroughly validated to simulate the acquisition of PA signals from complex media with a physical PA imaging system. The toolchain is based on two separate modeling steps: simulating the initial pressure field first and subsequently simulating the acoustic propagation of the generated PA waves. Furthermore, a new sampling strategy for the initial pressure was proposed to reduce discretization artifacts in the reconstructed images. The toolchain developed allows for quantitative simulations of PA signals, with a proper scaling factor for speckle intensities. The toolchain was designed to be as generic as possible, such that many different types of PA systems and media can be simulated.

The optical fluence was calculated by means of an MC method, which enables the simulation of the light distribution in complex media according to the RTE. A grid-based approach was taken to simulate the optical fluence, as this allowed for efficient raytracing of the photon packets and a convenient interpolation to the grid-based acoustic wave field solver of k-Wave. The developed MC software employs the microscopic Beer–Lambert law method for faster convergence and is parallelized for improved performance compared with traditional MC codes.²⁸^,³⁹ Although more sophisticated methods have been described in the literature that make use of graphical processing units to increase computational performance,⁴⁰ the developed software was mainly aimed to offer a practical and user-friendly modeling platform. A drawback of the used grid-based method, however, is a less accurate simulation of refraction and reflection effects, compared to mesh-based methods.⁴¹^–⁴³ The effect of refraction and reflection was deemed relatively small compared with scattering and absorption effects in most biological tissues.

The acoustic wave modeling software k-Wave was used in the toolchain to simulate the propagation of PA waves. k-Wave solves the wave equations governing the behavior of PA signal propagation through media and is therefore able to model all relevant acoustic phenomena, i.e., second and higher order scattering, frequency-dependent attenuation, reflection, and refraction. Due to the grid-based nature of k-Wave, staircasing occurs when discretizing the initial pressure field. Staircasing is a fundamental issue in grid-based discretization and has been studied in other physical domains as well, for example, in electromagnetics.⁴⁴ The effect of staircasing in band-limited pressure fields in PA was investigated in this study. The numerical method to prescribe $p_{0}$ in k-Wave was optimized to suppress staircasing errors in the simulated PA data. The preprocessing proposed step can achieve a similar performance as the direct simulation method (finer simulation grid) and the method introduced by Wise et al. (integration of point absorbers) in reducing the reconstructed intensity error, by decreasing the spatial spacing only when calculating p₀ and leaving k-Wave’s grid spacing intact. The staircasing artifacts at the disc’s boundary were absent in the Wise et al. method. Therefore, the Wise et al. method preserves the geometrical shape of the absorbers better, even at smaller sampling spacings. Still, it suffered from other artifacts in the region within the disc. These artifacts are probably caused by the BLI thresholding at 5% and by a slight nonuniform distribution of the integration points within the disc. Both the Wise et al. method and the direct simulation method led to a significant increase in computation time, as they were $˜ 17$ and 125 times slower than the preprocessing method, respectively, when simulating a 5-mm radius disc at a $p_{0}$ sampling of 40 PPWL. Note that the acoustic simulations were performed in 2D and that the required simulation time for the methods can differ in 3D. The preprocessing step therefore effectively improves the dynamic range of the PA images that can be simulated accurately. Furthermore, the preprocessing step is performed independent of the actual optical and acoustic simulations, allowing it to be used effectively with other methods and software tools than used in this study. To ensure a convenient workflow for the user, the optimal grid spacing in the preprocessing method can be straightforwardly determined using an empirical fit [Eq. (13)], after specifying the required dynamic range in the simulated PA images.

To validate the toolchain, simulated images of two phantoms and an ex vivo carotid plaque were compared with reconstructed images from a real PA set-up. The validation was performed on PA images, as this is a convenient and natural way to process and analyze the data. However, as the toolchain provides the raw channel data, the processing of the data is not limited to imaging only. Although all simulation results agree well with the real PA measurements in general, some small deviations are still present. The axial resolution of the simulated PSF was close to the in vitro measured resolution with only an absolute error of $40 μ m$ (0.2 wavelengths) as measured by the FWHM of the wire signal. The underestimation of FWHM of the wire in the simulation may be caused by the relatively large diameter of the wire ( $280 μ m$ ) which is larger than the wavelength at the probe’s center frequency ( $200 μ m$ ). Still, the simulated grating lobes agreed well with the experimental data in terms of reconstructed intensity and position. In this study, the system’s frequency response filter was defined by means of the probe’s center frequency and bandwidth using a Hann window in the temporal frequency domain. Although this approach is relatively straightforward and accurate for the in vitro system used, the PSF could be further improved by measuring the impulse response using a point absorber much smaller than the acoustic wavelength or directly using a calibrated hydrophone.

The simulated pressure amplitude and depth of maximum PA signal intensity in the reconstructed images matched well with the in vitro data for all concentrations of intralipid. The simulated pressure profiles as function of depth were accurate, although the pressure amplitude was generally underestimated in regions next to the peak. As this error is also present in the case of 0% intralipid, the error is probably caused by an inaccurate simulation of the light source. Furthermore, optical reflections in the in vitro set-up may also contribute to the discrepancy between the simulated and in vitro profiles. Moreover, it should be noted that, in the simulations, the Henyey–Greenstein phase function was used to calculate the optical scattering angle distribution, as is often done in tissue optics, whereas the scattering in intralipid is dominated by Mie scattering. These different phase functions may change the appearance of the pressure profiles slightly, even for equal reduced scattering coefficients and scattering anisotropy factors. Still the results suggest that simulations can replicate the actual optical wave propagation with high accuracy.

The PA speckle present in the ex vivo images was simulated using a Gaussian model to add random fluctuations to the initial pressure field. The appearance of PA speckle patterns in the simulated images, however, is highly dependent on the segmentation and prescription of the initial pressure. Although the Gaussian model gives reasonable results, other distributions or techniques to define the initial pressure may be investigated to increase the realism of the simulated PA signals in tissue models.

5. Conclusion

In this study, a flexible toolchain to simulate PA signals has been developed and validated. MC software was developed to determine the optical fluence in the illuminated medium and the open-source k-Wave toolbox was used to simulate the acoustic propagation of the pressure field over time. A computationally efficient preprocessing step was developed to decrease aliasing errors in the prescribed initial pressure field, which cause incorrect signals in the reconstructed PA images. It was shown that this preprocessing step could achieve a similar reduction in aliasing errors compared to the traditional approach of decreasing the spatial step in the k-Wave simulation. The toolchain was validated by in vitro measurements of a black wire, a three-channel phantom and ex vivo images of two human carotid plaques. It was shown that the PSF could be simulated accurately in turbid media. The overall appearance and SNR of the simulated PA images of the three-channel phantom and the two plaques showed good agreement to real PA measurements, indicating the great potential of the PA simulation in a more complex and heterogeneous media using the proposed toolchain.

6. Appendix

The energy of a 1D signal can be expressed using the following equation:

Eq. (14)

E = \sum_{n = 0}^{N - 1} s (n)^{2} Δ x,

where

s (n)

is the signal at position

n

, and

Δ x

is the spatial pixel spacing.

For a zero mean noise signal, the expected energy can also be written as function of the variance of $s$ :

Eq. (15)

E = N \cdot Var (s) \cdot Δ x .

To simulate a constant speckle intensity for all possible

Δ x

, the energy per unit length must be constant. Furthermore, because the total power of white noise scales proportionally with the bandwidth, the energy per unit of frequency should also be constant.

I_{speckle}

can be related to

E

by applying a proper scaling:

Eq. (16)

I_{speckle} = \frac{E}{L \cdot f_{s}} = \frac{E}{(Δ x \cdot N) \cdot \frac{1}{Δ x}} = \frac{E}{N} = Var (s) \cdot Δ x,

where

L

is the total length of the signal and

f_{s}

is the spatial sampling frequency. Finally,

I_{speckle}

can be related to the standard deviation of s,

σ_{s}

:

Eq. (17)

σ_{s} = \sqrt{Var (s)} = \sqrt{\frac{I_{speckle}}{Δ x}} .

The derivation can be generalized for multiple dimensions to obtain Eq. (6).

Disclosures

The authors declare no conflicts of interest.

Acknowledgments

We thank Camilo Cano from the Biomedical Engineering Department, Eindhoven University of Technology for his help with the acquisition of the experimental phantom and plaque data. This research was funded by the Horizon 2020 Framework Programme (731771), and Stichting Lijf en Leven (37).

Code, Data, and Materials Availability

The code and data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1.

T. Tarvainen et al., “Bayesian image reconstruction in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging, 32 (12), 2287 –2298 https://doi.org/10.1109/TMI.2013.2280281 ITMID4 0278-0062 (2013). Google Scholar

2.

S. Jeon and C. Kim, “Deep learning-based speed of sound aberration correction in photoacoustic images,” Proc. SPIE, 11240 112400J https://doi.org/10.1117/12.2543440 PSISDG 0277-786X (2020). Google Scholar

3.

J. Tick, A. Pulkkinen and T. Tarvainen, “Modelling of errors due to speed of sound variations in photoacoustic tomography using a Bayesian framework,” Biomed. Phys. Eng. Express, 6 (1), 015003 https://doi.org/10.1088/2057-1976/ab57d1 (2020). Google Scholar

4.

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 https://doi.org/10.1117/1.JBO.22.4.041008 JBOPFO 1083-3668 (2016). Google Scholar

5.

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 (4), 1468 https://doi.org/10.1364/BOE.7.001468 BOEICL 2156-7085 (2016). Google Scholar

6.

J. Gröhl et al., “Deep learning for biomedical photoacoustic imaging: a review,” Photoacoustics, 22 (January), 100241 https://doi.org/10.1016/j.pacs.2021.100241 (2021). Google Scholar

7.

S. Agrawal et al., “Modeling combined ultrasound and photoacoustic imaging: simulations aiding device development and artificial intelligence,” Photoacoustics, 24 100304 https://doi.org/10.1016/j.pacs.2021.100304 (2021). Google Scholar

8.

S. Manohar and D. Razansky, “Photoacoustics: a historical review,” Adv. Opt. Photonics, 8 (4), 586 https://doi.org/10.1364/AOP.8.000586 AOPAC7 1943-8206 (2016). Google Scholar

9.

L. V. Wang and H.-I. Wu, Biomedical Optics: Principles and Imaging, John Wiley & Sons( (2012). Google Scholar

10.

I. Krasnikov, A. Seteikin and B. Roth, “Advances in the simulation of light–tissue interactions in biomedical engineering,” Biomed. Eng. Lett., 9 (3), 327 –337 https://doi.org/10.1007/s13534-019-00123-x (2019). Google Scholar

11.

C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt., 18 (5), 050902 https://doi.org/10.1117/1.JBO.18.5.050902 JBOPFO 1083-3668 (2013). Google Scholar

12.

J. Gu and Y. Jing, “Modeling of wave propagation for medical ultrasound: a review,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 62 (11), 1979 –1992 https://doi.org/10.1109/TUFFC.2015.007034 ITUCER 0885-3010 (2015). Google Scholar

13.

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 https://doi.org/10.1117/1.3360308 JBOPFO 1083-3668 (2010). Google Scholar

14.

N. Akhlaghi et al., “Multidomain computational modeling of photoacoustic imaging: verification, validation, and image quality prediction,” J. Biomed. Opt., 24 (12), 121910 https://doi.org/10.1117/1.JBO.24.12.121910 JBOPFO 1083-3668 (2019). Google Scholar

15.

E. Aytac-Kipergil et al., “Development of a fiber laser with independently adjustable properties for optical resolution photoacoustic microscopy,” Sci. Rep., 6 38674 https://doi.org/10.1038/srep38674 SRCEC3 2045-2322 (2016). Google Scholar

16.

C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol., 54 (19), R59 –R97 https://doi.org/10.1088/0031-9155/54/19/R01 (2009). Google Scholar

17.

T. J. Allen and P. C. Beard, “Pulsed near-infrared laser diode excitation system for biomedical photoacoustic imaging,” Opt. Lett., 31 (23), 3462 https://doi.org/10.1364/OL.31.003462 OPLEDP 0146-9592 (2006). Google Scholar

18.

K. Daoudi et al., “Handheld probe integrating laser diode and ultrasound transducer array for ultrasound/photoacoustic dual modality imaging,” Opt. Express, 22 (21), 26365 https://doi.org/10.1364/OE.22.026365 OPEXFF 1094-4087 (2014). Google Scholar

19.

L. V. W. J. Yao, “Photoacoustic microscopy,” Laser Photonic Rev., 7 (5), 758 –778 https://doi.org/10.1002/lpor.201200060 (2013). Google Scholar

20.

E. Merčep, X. L. Deán-Ben and D. Razansky, “Imaging of blood flow and oxygen state with a multi-segment optoacoustic ultrasound array,” Photoacoustics, 10 48 –53 https://doi.org/10.1016/j.pacs.2018.04.002 (2018). Google Scholar

21.

A. Dima and V. Ntziachristos, “Non-invasive carotid imaging using optoacoustic tomography,” Opt. Express, 20 (22), 25044 https://doi.org/10.1364/OE.20.025044 OPEXFF 1094-4087 (2012). Google Scholar

22.

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

23.

G. Paltauf, P. R. Torke and N. Robert, “Modeling photoacoustic imaging with a scanning focused detector using Monte Carlo simulation of energy deposition,” J. Biomed. Opt., 23 (12), 121607 https://doi.org/10.1117/1.JBO.23.12.121607 JBOPFO 1083-3668 (2018). Google Scholar

24.

C. Fadden and S. R. Kothapalli, “A single simulation platform for hybrid photoacoustic and RF-acoustic computed tomography,” Appl. Sci., 8 (9), 1568 https://doi.org/10.3390/app8091568 (2018). Google Scholar

25.

A. Sharma et al., “Photoacoustic imaging depth comparison at 532-, 800-, and 1064-nm wavelengths: Monte Carlo simulation and experimental validation,” J. Biomed. Opt., 24 (12), 121904 https://doi.org/10.1117/1.JBO.24.12.121904 JBOPFO 1083-3668 (2019). Google Scholar

26.

J.-W. Muller et al., “Towards in vivo photoacoustic imaging of vulnerable plaques in the carotid artery,” Biomed. Opt. Express, 12 (7), 4207 https://doi.org/10.1364/BOE.430064 BOEICL 2156-7085 (2021). Google Scholar

27.

S. A. Prahl, “A Monte Carlo model of light propagation in tissue,” Proc. SPIE, 10305 1030509 https://doi.org/10.1117/12.2283590 PSISDG 0277-786X (1989). Google Scholar

28.

A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A, 29 (10), 2110 https://doi.org/10.1364/JOSAA.29.002110 JOAOD6 0740-3232 (2012). Google Scholar

29.

E. S. Wise et al., “Representing arbitrary acoustic source and sensor distributions in Fourier collocation methods,” J. Acoust. Soc. Am., 146 (1), 278 –288 https://doi.org/10.1121/1.5116132 JASMAN 0001-4966 (2019). Google Scholar

30.

M. U. Arabul et al., “Toward the detection of intraplaque hemorrhage in carotid artery lesions using photoacoustic imaging,” J. Biomed. Opt., 22 (4), 041010 https://doi.org/10.1117/1.JBO.22.4.041010 JBOPFO 1083-3668 (2016). Google Scholar

31.

R. Michels, F. Foschum and A. Kienle, “Optical properties of fat emulsions,” Opt. Express, 16 (8), 5907 https://doi.org/10.1364/OE.16.005907 OPEXFF 1094-4087 (2008). Google Scholar

32.

H. Buiteveld, J. H. M. Hakvoort and M. Donze, “Optical properties of pure water,” Proc. SPIE, 2258 174 –183 https://doi.org/10.1117/12.190060 PSISDG 0277-786X (1994). Google Scholar

33.

M. U. Arabul et al., “Investigation on the effect of spatial compounding on photoacoustic images of carotid plaques in the in vivo available rotational range,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 65 (3), 440 –447 https://doi.org/10.1109/TUFFC.2018.2792903 ITUCER 0885-3010 (2018). Google Scholar

34.

L. Peralta et al., “Coherent multi-transducer ultrasound imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 66 (8), 1316 –1330 https://doi.org/10.1109/TUFFC.2019.2921103 (2019). Google Scholar

35.

M. J. C. van Gemert et al., “Optical properties of human blood vessel wall and plaque,” Lasers Surg. Med., 5 (3), 235 –237 https://doi.org/10.1002/lsm.1900050305 LSMEDI 0196-8092 (1985). Google Scholar

36.

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

37.

M. D. Verweij et al., “Simulation of ultrasound fields,” Comprehensive Biomedical Physics, 2 Elsevier( (2014). Google Scholar

38.

V. Periyasamy and M. Pramanik, “Advances in Monte Carlo simulation for light propagation in tissue,” IEEE Rev. Biomed. Eng, 10 122 –135 https://doi.org/10.1109/RBME.2017.2739801 (2017). Google Scholar

39.

S. L. Jacques and L. Wang, “Monte Carlo Modeling of light transport in tissues,” Opt. Response Laser-Irradiated Tissue, 2607 (713), 73 –100 https://doi.org/10.1007/978-1-4757-6092-7_4 (1995). Google Scholar

40.

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 https://doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 (2009). Google Scholar

41.

Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express, 1 (1), 165 –175 https://doi.org/10.1364/BOE.1.000165 BOEICL 2156-7085 (2010). Google Scholar

42.

A. A. Leino, A. Pulkkinen and T. Tarvainen, “ValoMC: a Monte Carlo software and MATLAB toolbox for simulating light transport in biological tissue,” OSA Contin., 2 (3), 957 https://doi.org/10.1364/OSAC.2.000957 (2019). Google Scholar

43.

S. Yan, A. P. Tran and Q. Fang, “Dual-grid mesh-based Monte Carlo algorithm for efficient photon transport simulations in complex three-dimensional media,” J. Biomed. Opt., 24 (2), 020503 https://doi.org/10.1117/1.JBO.24.2.020503 JBOPFO 1083-3668 (2019). Google Scholar

44.

J. A. McCammon, “Dielectric boundary smoothing in finite difference solutions of the Poisson equation,” J. Comput. Chem., 12 (7), 909 –912 https://doi.org/10.1002/jcc.540120718 JCCHDD 0192-8651 (1991). Google Scholar

Biography

Jan-Willem Muller received his BSc and MSc degrees in biomedical engineering from TU/e, in 2016 and 2018, respectively. He is currently pursuing his PhD at the Photoacoustic and Ultrasound Laboratory Eindhoven. Currently, his research focuses on photoacoustic and ultrasound imaging of atherosclerotic plaques in the carotid artery. His research interests include in silico modeling of photoacoustic and ultrasound imaging, beamforming methods, strain imaging, and displacement estimate regularization.

Mustafa Ü. Arabul received his BSc degree in electrical and electronics engineering from the Middle East Technical University, Ankara, Turkey, in 2011. Next, he joined the Medical Imaging Lab at the Institute of Biomedical Engineering, Bogaziçi University, Istanbul, and achieved his MSc degree in biomedical engineering in 2013. He pursued his PhD in the Photoacoustics and Ultrasound Laboratory of Eindhoven (PULS/e Lab), in the Cardiovascular Biomechanics group of the Department of Biomedical Engineering at Eindhoven University of Technology, The Netherlands. His research focused on photoacoustic imaging of carotid arteries and fundamental characterization and preclinical validation of photoacoustic imaging (2018). Currently, he is working as a research scientist at ASML with the focus on semiconductor metrology.

Hans-Martin Schwab received his bachelor’s degree in engineering science from the TU Berlin, in 2011. In 2013, he finished his master’s degree in electrical engineering with a focus on medical imaging at Ruhr University Bochum, where he also obtained his PhD on photoacoustic imaging in 2018. After this, he became a postdoctoral researcher at the PULS/e at Eindhoven University of Technology, where he was appointed as an assistant professor in 2021. His current interests focus on approaches for ultrasound acquisition and model-based image reconstruction.

Marcel C. M. Rutten is an assistant professor at the TU/e Department of Biomedical Engineering, Research Group Cardiovascular Biomechanics. His research focuses on vulnerable plaques, their characterization by means of techniques, such as ultrasound and photoacoustics, and the development of biomechanical models.

Marc R. H. M. van Sambeek is a part-time professor of application of engineering in vascular surgery and endovascular therapy at the Cardiovascular Biomechanics Research Group. His main affiliation is with the Catharina Hospital Eindhoven, as a vascular surgeon. He specializes in aortic aneurysm, carotid artery stenosis, and thoracic outlet syndrome, and he is interested in minimal invasive vascular surgery. At TU/e, his research focuses on models using patient-specific imaging to improve clinical decision making and intervention.

Min Wu is an assistant professor at PULS/e Lab at TU/e. Her research area focuses on photoacoustic imaging. Her research interests include the development of photoacoustic imaging system and characterization of different tissue compositions with multispectral photoacoustic imaging. Her research aims to support clinical diagnosis and treatment in cardiovascular and orthopedic applications.

Richard G. P. Lopata received his MSc degree in biomedical engineering in 2004 (TU/e) and his PhD from the Radboudumc, on “2-D and 3-D ultrasound strain imaging: methods and in vivo applications.” He has been an associate professor (TU/e), heading the PULS/e Lab research group since 2014, which facilitates research in the areas of ultrasound functional imaging, photoacoustics, and image-based modeling aimed to facilitate/improve clinical decision making for cardiovascular, musculoskeletal, and abdominal applications.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International 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

Jan-Willem Muller, Mustafa Umit Arabul, Hans-Martin Schwab, Marcel C. M. Rutten, Marc R. H. M. van Sambeek, Min Wu, and Richard G. P. Lopata "Modeling toolchain for realistic simulation of photoacoustic data acquisition," Journal of Biomedical Optics 27(9), 096005 (14 September 2022). https://doi.org/10.1117/1.JBO.27.9.096005

Received: 7 June 2022; Accepted: 31 August 2022; Published: 14 September 2022

Access the abstract

JOURNAL ARTICLE
17 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 7 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Acoustics

Data modeling

In vitro testing

Speckle

Tissue optics

Photoacoustic spectroscopy

Tissues

1.

Introduction

2.

Methods

2.1.

Optical Fluence Modeling

Eq. (1)

Eq. (2)

Eq. (3)

2.2.

Acoustic Simulation and Sensor Definition

2.3.

Initial Pressure Definition

Eq. (4)

Fig. 1

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Fig. 2

2.4.

Staircasing Error Reduction Analysis

Eq. (10)

Eq. (11)

Eq. (12)

2.5.

In Vitro and In Silico Experimental Setup

Fig. 3

2.6.

Verification of Simulated PSF in Turbid Media

2.7.

Three-Channel Phantom

Fig. 4

2.8.

Ex Vivo Plaque Imaging

Table 1

3.

Results

3.1.

Initial Pressure Preprocessing

Fig. 5

Eq. (13)

Fig. 6

3.2.

Fluence Simulation in Turbid Medium

Fig. 7

3.3.

Three-Channel Phantom

Fig. 8

3.4.

Ex Vivo Plaque Imaging

Fig. 9

4.

Discussion

5.

Conclusion

6.

Appendix

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Disclosures

Acknowledgments

Code, Data, and Materials Availability

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years