Linking information theory and thermodynamics to spatial resolution in photothermal and photoacoustic imaging

Burgholzer, P.; Mayr, G.; Thummerer, G.; Haltmeier, M.

doi:10.1063/5.0023986

In this Tutorial, we combine the different scientific fields of information theory, thermodynamics, regularization theory, and non-destructive imaging, especially for photoacoustic and photothermal imaging. The goal is to get a better understanding of how information gaining for subsurface imaging works and how the spatial resolution limit can be overcome by using additional information. Here, the resolution limit in photoacoustic and photothermal imaging is derived from the irreversibility of attenuation of the pressure wave and of heat diffusion during the propagation of the signals from the imaged subsurface structures to the sample surface, respectively. The acoustic or temperature signals are converted into so-called virtual waves, which are their reversible counterparts and which can be used for image reconstruction by well-known ultrasound reconstruction methods. The conversion into virtual waves is an ill-posed inverse problem, which needs regularization. The reason for that is the information loss during signal propagation to the sample surface, which turns out to be equal to the entropy production. As the entropy production from acoustic attenuation is usually small compared to the entropy production from heat diffusion, the spatial resolution in acoustic imaging is higher than in thermal imaging. Therefore, it is especially necessary to overcome this resolution limit for thermographic imaging by using additional information. Incorporating sparsity and non-negativity in iterative regularization methods gives a significant resolution enhancement, which was experimentally demonstrated by one-dimensional imaging of thin layers with varying depth or by three-dimensional imaging, either from a single detection plane or from three perpendicular detection planes on the surface of a sample cube.

I. INTRODUCTION

Metrology can be described as collecting information from samples. Imaging techniques provide much information, as is indicated by the saying “a picture is worth a thousand words.” Non-destructive evaluation (NDE) or biomedical imaging often image deep structures in the samples' interior without destroying the samples. Therefore, an information theoretical viewpoint is instructive to determine the amount of information on subsurface or other embedded structures, which can be gained by measurements on the sample surface. Information theory is a branch of applied mathematics, electrical engineering, and computer science involving the quantification of information. Information theory was developed by Shannon¹ to find fundamental limits on signal processing operations such as compressing data and on reliably storing and communicating data. Information processing is a physical activity, which has to obey the laws of (non-equilibrium) thermodynamics. This was originally recognized by Szilard,² and Landauer³ captured it with his aphorism: “Information is physical.”

Imaging of subsurface features from data measured on the sample surface usually results in an ill-posed or an ill-conditioned inverse problem, which needs regularization,^4–6 as illustrated in Fig. 1 for photoacoustic and photothermal imaging. Regularization typically involves additional assumptions, such as the smoothness of the solution. Usually, scientists working on the inverse problem of subsurface imaging are barely working on non-equilibrium thermodynamics and information theory and vice versa. This Tutorial should enable researchers from both sides to get more insight from the counterpart. The models used for subsurface imaging and for thermodynamics are kept simple in the beginning, allowing to investigate the essential connections between thermodynamics and regularization of inverse problems. As a first step, a one-dimensional (1D) model without boundaries and a Dirac delta-like structure at a varying depth is used for imaging, either inducing thermal waves described by heat diffusion or acoustic pressure signals showing acoustic attenuation as in a liquid. For describing the relevant results from non-equilibrium thermodynamics, we follow a paper from Esposito and van den Broeck about the “Second law and Landauer principle far from equilibrium,”⁷ where we describe the photothermal or photoacoustic measurement process in terms of thermodynamics. The description can be simplified as the Hamiltonian does not change in time. A short laser light pulse (Dirac delta-like) is used for excitation, which brings the sample suddenly out of thermal equilibrium by heating the structures via optical absorption and thereby generating an initial pressure distribution, and then the system returns slowly to equilibrium either by heat diffusion (“thermal”) or by attenuation of the pressure wave (“acoustic”).

FIG. 1.

View large Download slide

Schematic sketch of the imaging of absorbers within a sample inspired by Li et al.⁸ (a) Sample containing absorbers as an internal structure, which should be imaged by suddenly heating this structure, e.g., by absorbing scattered light from a laser pulse. (b) Propagation of induced signals, such as thermal or pressure waves to the sample surface, where they are detected. Signals are blurred due to heat diffusion or acoustic attenuation, respectively. (c) This signal blurring, which reflects noise and entropy production (Sec. III), can be compensated to a certain extent when calculating the virtual wave. The virtual wave is a solution of the ideal wave equation without diffusion or attenuation. Calculating the virtual wave is an ill-posed inverse problem that needs regularization [Eq. (1)]. (d) Image reconstruction by wave backpropagation, either of the virtual waves gained from the temperature signals in Photo-Thermal Computed Tomography (PTCT) or of the attenuation compensated virtual acoustic waves in Photo-Acoustic Computed Tomography (PACT).

After presenting these simple models for the wave propagation and the thermodynamic approach, they are linked in the frequency domain to evaluate the physical background of regularization. Later on, more realistic samples in higher dimensions and with realistic boundary conditions are described. But before stepping into the detail, in the remaining part of Sec. I, previous work in subsurface imaging and in thermodynamics is described.

In former times, images always have shown what could be detected with human eyes. Looking under the sample surface in a non-destructive way was possible only for optically transparent samples. This situation changed drastically with the discovery of x rays by Röntgen in 1895. Nowadays, the interaction of sample structures with electromagnetic waves in the whole frequency range from eddy current (EC), radar, terahertz (THz), infrared (IR), up to UV, and x-ray radiation is used for imaging, but also with particles or elastic (acoustic) waves. The amount of information, which can be gained by measurements on the sample surface, depends on the interaction of these waves or particles with structures of the sample, such as scattering or absorption, and on the detector (e.g., sensitivity or bandwidth). For interior structures, the available information for imaging is limited also by the information loss during wave propagation from the imaged structures to the sample’s surface, caused by diffusion, dissipation, or scattering (Fig. 1 or Burgholzer and Hendorfer⁹). A unifying framework for treating diverse diffusion-related periodic phenomena under the global mathematical label of diffusion-wave fields has been developed by Mandelis,¹⁰ including thermal waves, charge-carrier-density waves, and diffuse-photon-density waves, but also modulated eddy currents, neutron waves, or harmonic mass-transport diffusion waves.

There have been made several attempts to compensate the diffusion, dissipation, or scattering during wave propagation to get a higher resolution for the reconstructed images of the samples interior. We have shown that thermodynamical fluctuations are the cause for an entropy production, which is equal to the information loss and limits this compensation.⁹ Therefore, also the spatial resolution for non-destructive imaging (NDI) at a certain depth is limited. For waves that satisfy the wave equation, no information is lost during the propagation, because the time reversed wave is a solution of the same wave equation. The inverse problem of “backprojection” can be exactly solved by, e.g., time reversal.¹¹ It is called a virtual wave signal, because in reality always dissipation, diffusion, or scattering causes a loss of information during propagation. For attenuated acoustic waves^12–24 and thermal diffusion,^9,25–30 a linear relation between the time discretized measured signal $Y_{meas}$ and the virtual wave signal $Y_{virt}$ ⁠, which is a solution of the wave equation, was established and can be written as

Y_{meas} = M Y_{virt} .

(1)

The matrix $M$ depends on the kind of wave propagation and was calculated already for acoustic waves having a power law attenuation and for heat diffusion. The solution of Eq. (1) constitutes an ill-posed or ill-conditioned inverse problem. The direct inversion of this equation would, therefore, cause severe noise amplification resulting in useless solutions. Therefore, regularization methods are used for the inversion. This is a consequence of the information loss during signal propagation to the surface. $Y_{meas}$ could be the measured acoustic pressure or temperature, and $Y_{virt}$ could be the virtual pressure or temperature, respectively.

Equation (1) was derived for attenuated acoustic waves and thermal diffusion, but the concept of deriving the resolution limit from entropy production seems to be useful for non-destructive and biomedical imaging, in general. For example, applications of a transformation of the diffusive electromagnetic wave into a wave field were shown by Lee et al.^31,32 and Gershenson³³ for geophysical inverse problems. The physical reason for the entropy production is fluctuations, which is the noise added to an ideally noise-free signal and can be statistically described by stochastic processes.⁹ The ratio between the ideal signal amplitude and the noise amplitude is called the signal-to-noise ratio (SNR) and plays an essential role in this Tutorial. Noise is identified as the thermodynamic fluctuations around a certain time varying mean value and is shown in Sec. III to be the mechanism for entropy production and information loss.

One comprehensive letter about the relevant non-equilibrium thermodynamics, the second law, and the connection between entropy production and information loss was published in 2011 by Esposito and van den Broeck.⁷ They showed that for two different non-equilibrium states evolving to the same equilibrium state, the entropy production $Δ_{i} S$ during the evolution from one state to the other is equal to the information loss $Δ I = k_{B} Δ D$ ⁠, where $k_{B}$ is the Boltzmann constant and $Δ D$ is the difference of the Kullback–Leibler divergence $D$ ⁠, also called relative entropy of these states. $D$ is a measure of how “far” a certain state is away from equilibrium.³⁴ The entropy production for the macroscopic states with small fluctuations around equilibrium turns out to be, in a good approximation, equal to the dissipated energy $Δ Q$ ⁠, which is the dissipated heat or the heat transported by diffusion, divided by the mean temperature $T_{mean}$ ⁠, so $Δ_{i} S = Δ Q / T_{mean} = k_{B} Δ D$ ⁠.²⁶

The thermodynamic fluctuations in macroscopic samples are usually so small that they can be neglected—but not for ill-posed inverse problems such as described by the inversion of Eq. (1). The fluctuations are highly amplified due to the ill-posed problem of image reconstruction. As long as the macroscopic mean-value equations describe the mean entropy production, the resolution limit depends only on the amplitude of the fluctuations and not on the actual stochastic process including all correlations.²⁶ Therefore, we have used in the past a simple Gauss–Markov process to describe the measured signal as a time-dependent random variable.⁹ In addition, we have performed a preliminary work to describe heat diffusion as a Wiener process and have compared for a simulated signal the spatial resolution limits derived in temporal and spatial frequency domains, $ω$ - and $k$ -space, respectively.³⁵ Experimentally, these theoretical resolution limits were verified by thermographic reconstructions of inclined steel rods in an epoxy sample, heated by a short light pulse²⁷ (Sec. IV B 1).

II. PHOTOACOUSTIC AND PHOTOTHERMAL IMAGING MODEL

In photothermal and photoacoustic imaging, a short light pulse, e.g., from a laser, is used to heat subsurface structures by the absorbed (scattered) light, such as blood vessels in tissue or light absorbing structures in epoxy resin, creating a temperature distribution $T_{0} (r)$ at location $r$ immediately after the pulse [Fig. 1(a)]. This sudden temperature increase causes an initial pressure distribution $p_{0} (r)$ proportional to the absorbed optical energy density $A (r) = T_{0} (r) C_{p} ρ_{0}$ ⁠, with a dimensionless material constant $Γ = β c^{2} / C_{p}$ ⁠, the Grüneisen coefficient.^36–38 The initial pressure distribution can be written as

p_{0} (r) = Γ A (r) = \frac{β c^{2}}{C_{p}} A (r) = β c^{2} ρ_{0} T_{0} (r),

(2)

where $C_{p}$ is the specific heat capacity, $ρ_{0}$ is the ambient density, $c$ is the sound speed, and $β$ is the thermal volume expansion coefficient at constant pressure. The initial pressure distribution causes a virtual pressure wave $p_{virt} (r, t)$ as a function of space $r$ and time $t$ ⁠. This wave is called “ideal” as no acoustic attenuation or dispersion during propagation is assumed and is “virtual” as attenuation and dispersion are always present for real waves, even if acoustic attenuation can often be neglected for lower frequencies. The wave equation describes the acoustic pressure $p_{virt} (r, t)$ as a function of space $r$ and time $t$ and can be written as^37,38

(\nabla^{2} - \frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}}) p_{virt} (r, t) = - \frac{1}{c^{2}} \frac{\partial}{\partial t} p_{0} (r) δ (t),

(3)

where $\nabla^{2}$ is the Laplacian (second derivative in space). The source term on the right side of Eq. (3) ensures that the pressure just after the short excitation pulse, which is modeled by the temporal Dirac delta function $δ (t)$ ⁠, is the initial pressure distribution $p_{0} (r)$ ⁠. The same equation is valid for a virtual temperature wave, which is defined as $T_{virt} := p_{virt} (r, t) / (β c^{2} ρ_{0})$ and gives

(\nabla^{2} - \frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}}) T_{virt} (r, t) = - \frac{1}{c^{2}} \frac{\partial}{\partial t} T_{0} (r) δ (t) .

(4)

The virtual wave is the induced virtual pressure wave, but multiplied by a material constant to get a temperature measured in Kelvin. As it has no direct physical representation but is a mathematical model for thermographic reconstruction, it is called “virtual temperature wave.”

The actual temperature evolution $T (r, t)$ can be described by the heat diffusion equation³⁹

(\nabla^{2} - \frac{1}{α} \frac{\partial}{\partial t}) T (r, t) = - \frac{1}{α} T_{0} (r) δ (t),

(5)

where $T (r, t)$ is the temperature as a function of space and time, and $α$ is the thermal diffusivity, which is assumed to be homogeneous in the sample. This description of the thermal diffusion is based on Fourier’s law, which is valid for macroscopic samples, where the propagation distance is much larger than the phonon mean free path.⁴⁰

We denote by $\tilde{T} (r, ω)$ the temperature signal in the frequency domain, the $ω$ -space, which is related to $T (r, t)$ by the temporal Fourier transform

\tilde{T} (r, ω) = \int_{- \infty}^{\infty} T (r, t) \exp (i ω t) d t,

(6a)

and its inverse

T (r, t) = \frac{1}{2 π} \int_{- \infty}^{\infty} \tilde{T} (r, ω) \exp (- i ω t) d ω .

(6b)

Similarly, ${\tilde{p}}_{virt} (r, ω)$ and ${\tilde{T}}_{virt} (r, ω)$ denote the Fourier transforms of the virtual waves $p_{virt} (r, t)$ and $T_{virt} (r, t)$ ⁠, respectively. Taking the Fourier transform according to Eq. (6) of Eqs. (3)–(5) and using $δ (t) = \frac{1}{2 π} \int_{- \infty}^{\infty} \exp (- i ω t) d ω$ results in the Helmholtz equations

(\nabla^{2} + k {(ω)}^{2}) {\tilde{p}}_{virt} (r, ω) = \frac{i ω}{c^{2}} p_{0} (r),

(7)

(\nabla^{2} + k {(ω)}^{2}) {\tilde{T}}_{virt} (r, ω) = \frac{i ω}{c^{2}} T_{0} (r),

(8)

(\nabla^{2} + σ {(ω)}^{2}) \tilde{T} (r, ω) = - \frac{1}{α} T_{0} (r),

(9)

with $k (ω) \equiv ω / c$ and $σ {(ω)}^{2} \equiv i ω / α$ ⁠. For the virtual pressure wave and the virtual temperature wave, the wavenumber $k (ω)$ is real; for the temperature, the wavenumber $σ (ω)$ is complex with a real and imaginary part of equal size to describe the diffusive effect.

Acoustic attenuation and dispersion can also be described by a complex wavenumber, but usually the imaginary part is much smaller than the real part. As Stokes⁴¹ could already show in 1845 and later on, e.g., Shutilov,⁴² fundamental responses of the material system cause a relaxation time between pressure and density changes, which result in an attenuation of a propagating acoustic wave. For liquids, the acoustical absorption $α (ω)$ increases with the square of the angular frequency $ω$ ⁠, that is $α (ω) = α_{0} ω^{2}$ ⁠, with the material constant $α_{0}$ taking into account, e.g., the viscosity of the liquid. The pressure caused by a monochromatic acoustic plane wave of frequency $ω$ in a uniform attenuating medium at a point $x$ and at an instant $t$ can be expressed as

\begin{aligned} p_{ω} (x, t) & = \exp (i (k (ω) x - ω t) \exp (- α_{0} ω^{2} x) \\ \equiv \exp (i (K (ω) x - ω t)), \end{aligned}

(10)

where the acoustic attenuation coefficient $α (ω)$ can be written as the imaginary part of a complex wavenumber $K (ω) = k (ω) + i α (ω) = ω / c + i α_{0} ω^{2}$ ⁠. For the attenuated acoustic wave, the Helmholtz equation

(\nabla^{2} + K {(ω)}^{2}) \tilde{p} (r, ω) = \frac{i ω}{c^{2}} p_{0} (r)

(11)

is the same as for the virtual pressure wave in Eq. (7), but with the complex wavenumber.^13–15 The source term on the right side of Eq. (11) is the same as for the virtual pressure wave ${\tilde{p}}_{virt} (r, ω)$ in Eq. (7) and therefore a direct relation can be derived between the attenuated pressure wave and the virtual pressure wave on the same location $r$ ⁠.^13–15 Replacing $ω$ by $c K (ω)$ in Eq. (7) results in

(\nabla^{2} + K {(ω)}^{2}) {\tilde{p}}_{virt} (r, c K (ω)) = \frac{i c K (ω)}{c^{2}} p_{0} (r),

(12)

which is equal to Eq. (11), if multiplied by $ω / c K (ω)$ ⁠, and

\tilde{p} (r, ω) = \frac{ω}{c K (ω)} {\tilde{p}}_{virt} (r, c K (ω)) .

(13)

This is the sought relation between the attenuated pressure signal and the virtual pressure signal. The location $r$ is the same for ${\tilde{p}}_{virt}$ and $\tilde{p}$ ⁠, and the relation is valid in all dimensions. Transformation back from $ω$ –space to the time domain with the inverse Fourier transformation [Eq. (6b)] results in

p (r, t) = \frac{1}{2 π} \int_{- \infty}^{\infty} \frac{ω}{c K (ω)} {\tilde{p}}_{virt} (r, c K (ω)) \exp (- i ω t) d ω,

with

{\tilde{p}}_{virt} (r, c K (ω)) = \int_{- \infty}^{\infty} p_{virt} (r, t^{'}) \exp (i c K (ω) t^{'}) d t^{'},

(14)

which can be written as

\begin{aligned} p (r, t) & = \int_{- \infty}^{\infty} p_{virt} (r, t^{'}) M_{p} (t, t^{'}) d t^{'}, \\ with \\ M_{p} (t, t^{'}) & \equiv \frac{1}{2 π} \int_{- \infty}^{\infty} \frac{ω}{c K (ω)} \exp (i c K (ω) t^{'}) \exp (- i ω t) d ω \\ \approx \frac{1}{2 π} \int_{- \infty}^{\infty} \exp (i c K (ω) t^{'}) \exp (- i ω t) d ω \\ = \frac{1}{2 π} \sqrt{\frac{π}{α_{0} c t^{'}}} \exp (- \frac{{(t - t^{'})}^{2}}{4 α_{0} c t^{'}}) for t^{'} > 0. \end{aligned}

(15)

The factor $ω / (c K (ω))$ turns out to be approximately one for relevant frequencies and attenuation coefficients. This is shown for glycerine, a liquid with rather high viscosity, in Subsection II A: even for a high frequency of 100 MHz, the factor is 0.994. For lower frequencies or lower acoustic attenuation than in glycerine, this factor is even closer to one. Therefore, $M_{p} (t, t^{'})$ is in a good approximation to a Gaussian pulse, where the amplitude is reduced and the width is increased by a factor of $\sqrt{α_{0} c t^{'}}$ ⁠. For a thin absorbing planar layer in an infinite medium, which can be modeled as a Dirac delta function in the thickness coordinate $z$ ⁠, the virtual pressure $p_{virt} (z, t)$ for $z > 0$ is proportional to $δ (t - z / c)$ ⁠. From Eq. (15), one gets for the attenuated pressure wave

\begin{aligned} p (z, t) & \propto \int_{- \infty}^{\infty} δ (t^{'} - z / c) M_{p} (t, t^{'}) d t^{'} \\ \approx \frac{1}{2 π} \sqrt{\frac{π}{α_{0} z}} \exp (- \frac{{(t - z / c)}^{2}}{4 α_{0} z}), \end{aligned}

(16)

which is the inverse Fourier transform [Eq. (6b)] of $\exp (i K (ω) z)$ ⁠.⁴³ This is also true for a general complex $K (ω)$ ⁠, e.g., with different power laws than the square of the frequency, such as for porcine fat tissue with a power law exponent of 1.5 (Sec. IV A 1). More details on compensation of acoustic attenuation can be found in Refs. 12–24.

Equation (15) can be discretized to produce a matrix equation

p = M_{p} p_{virt},

(17)

where $p$ and $p_{virt}$ are the vectors of the attenuated and virtual pressure signal at discrete time steps, respectively. $M_{p}$ is the matrix at these time steps calculated from Eq. (15).

In a similar way, we have calculated a local relation between the actual temperature $T (r, t)$ and the virtual wave signal $T_{virt} (r, t)$ from Eqs. (4) and (5), or in the frequency domain from the Helmholtz equations (8) and (9).²⁵ Replacing $ω$ by $c σ (ω)$ in Eq. (8) gives Eq. (9) by multiplying with $\frac{i c}{α σ (ω)}$ and identifying

\tilde{T} (r, ω) = \frac{i c}{α σ (ω)} {\tilde{T}}_{virt} (r, c σ (ω)) .

(18)

The inverse Fourier transformation in Eq. (6b) can be integrated analytically and can be written as²⁵

\begin{aligned} T (r, t) & = \int_{- \infty}^{\infty} T_{virt} (r, t^{'}) M_{T} (t, t^{'}) d t^{'}, \\ with \\ M_{T} (t, t^{'}) & \equiv \frac{c}{\sqrt{π α t}} \exp (- \frac{c^{2} t^{' 2}}{4 α t}) for t > 0. \end{aligned}

(19)

Equation (19) like Eq. (18) connects the temperature signal to the virtual wave signal at the same location $r$ ⁠, but in the time domain instead of the temporal frequency domain.

For the example of a thin absorbing layer in an infinite medium, the virtual temperature wave consists of two Dirac delta pulses running in the positive and negative $z$ -direction proportional to $\frac{1}{2} δ (z \pm c t) = \frac{1}{2 c} δ (t \pm \frac{z}{c})$ ⁠. From Eq. (19), we get the 1D Green’s function

\begin{aligned} T (z, t) & \propto \int_{- \infty}^{\infty} \frac{1}{2} δ (z \pm c t^{'}) M_{T} (t, t^{'}) d t^{'} \\ = \frac{1}{\sqrt{4 π α t}} \exp (- \frac{z^{2}}{4 α t}) . \end{aligned}

(20)

Equation (19) can be discretized to produce a matrix equation

T = M_{T} T_{virt},

(21)

similar to Eq. (17), Eqs. (21) and (17) can be inverted only with appropriate regularization, as the matrices $M_{T}$ and $M_{p}$ are rank deficient. The truncated-singular value decomposition (T-SVD) method is used in the following to get the reconstructed signal $T_{rec}$ as an estimate for $T_{virt}$ from the thermographic signal $T$ or $p_{rec}$ as an estimate for $p_{virt}$ from the pressure signal $p$ ⁠, respectively. The truncation value for the smallest singular values is $1 / SNR$ ⁠, with SNR being the signal-to-noise ratio for the measured temperature or pressure.

Image reconstruction is now a two-stage process, where two successive inverse problems are solved. In photoacoustic imaging, in a first step, $p_{virt}$ from the measured pressure signal $p$ on the sample surface is determined and then in a second step the initial pressure distribution $p_{0} (r)$ is determined by acoustic reconstruction methods, such as backpropagation or time reversal reconstruction algorithms.¹¹ In photothermal imaging, the initial temperature distribution $T_{0} (r)$ is determined by the same reconstruction methods from the virtual temperature wave $T_{virt}$ ⁠, which was calculated in the first step from the measured surface temperature $T$ ⁠. This two-stage process for imaging can be used in 1D, 2D, and 3D. In one dimension, the second step, the reconstruction by backprojection, is rather trivial: the time signal of the virtual pressure or temperature wave multiplied by the constant sound speed $c$ gives directly the initial pressure or temperature at the depth $z = c t$ ⁠. Therefore, in the beginning, we will show a 1D example.

A. Compensation of acoustic attenuation for 1D photoacoustic imaging

In the frequency domain, according to Eq. (10), the amplitude of the wave component with frequency $ω$ is damped by the factor $\exp (- α_{0} ω^{2} z)$ during propagation from depth $z$ to the sample surface. For frequencies larger than the truncation frequency $ω_{cut}$ ⁠, the amplitude of these wave components is damped below the noise level, and at $ω_{cut}$ ⁠, we get

SNR \exp (- α_{0} ω_{cut}^{2} z) = 1 or ω_{cut} = \sqrt{\frac{\ln (SNR)}{α_{0} z}} .

(22)

For the spatial resolution in photoacoustic imaging, the minimal possible width of the acoustic signal in the time domain is essential. A small width enables high spatial resolution, which corresponds to a broad frequency bandwidth. If the frequency bandwidth is limited according to Eq. (22), the spatial resolution limit according to Nyquist is half the wavelength at this frequency

δ_{resolution} = \frac{π}{ω_{cut}} c = π c \sqrt{\frac{α_{0} z}{\ln (SNR)}} .

(23)

The resolution limit can be validated by the reconstruction of a thin absorbing planar layer in an infinite medium, which is modeled as a Dirac delta function in the thickness coordinate $z$ ⁠, and where the virtual pressure $p_{virt} (z, t)$ for $z > 0$ is proportional to $δ (t - z / c)$ ⁠. The Fourier transformation of a delta pulse is $\exp (i ω z / c)$ and shows all frequencies with equal amplitude. The reconstruction in the attenuation limited frequency bandwidth is the inverse Fourier transformation, where the frequency integral is taken from $- ω_{cut}$ to $+ ω_{cut}$ ⁠, which gives

\begin{aligned} p_{rec} (z, t) & = \frac{1}{2 π} \int_{- ω_{cut}}^{+ ω_{cut}} \exp (i ω \frac{z}{c}) \exp (- i ω t) d ω \\ = \frac{1}{π} \frac{1}{z / c - t} \sin (ω_{cut} (\frac{z}{c} - t)) . \end{aligned}

(24)

This is a sinc function, where the maximum at a distance $z$ is at $t = z / c$ ⁠, which is the arrival time of the virtual wave. The zero points are at $t = \frac{z}{c} \pm \frac{π}{ω_{cut}} = \frac{z \pm δ_{resolution}}{c}$ with $δ_{resolution}$ from Eq. (23). Compared to the measured pressure pulse without any compensation of attenuation given in Eq. (16), the width of the reconstructed pressure pulse, Eq. (24), is reduced by a factor of $\sqrt{\ln (SNR)}$ ⁠. A higher SNR allows a better resolution for the reconstruction.

Glycerine is a liquid with rather high viscosity. Acoustic attenuation in glycerine is about 100 times higher than in water. At a temperature of $25^{°} C$ and a frequency of 1 MHz, the sound velocity is $c = 1923 m / s$ and the attenuation is $2.5 m^{- 1}$ ⁠, which gives 22 dB/m or 0.22 dB/cm [multiplying the attenuation with $20 \log_{10}$ (⁠ $e$ ⁠)].⁴² This gives $α_{0} = α / ω^{2} = 2.5 / {(2 π 10^{6})}^{2} s^{2} / m = 6.3 \times 10^{- 14} s^{2} / m$ ⁠. Even for high acoustic attenuation as in glycerine and for high frequencies, e.g., up to 100 MHz, $ω / (c K (ω))$ can still be approximated by one in Eq. (16), because $c α_{0} ω = 0.076$ is small compared to one.

We will start with a simple 1D imaging problem: a thin planar layer at a certain depth $z$ below the surface of our detection plane. Due to acoustic attenuation, the rectangular pressure pulse gets a different shape according to Eq. (15), calculated by the discretized version of the matrix equation in (17). In our example, we take three layers at a depth $z$ of 1 mm, 5 mm, and 10 mm, with a layer thickness of 0.2 mm. Table I shows the truncation frequencies according to Eq. (22) and the resulting resolution limits from Eq. (23). Figure 2 shows the initial pressure distribution, which corresponds in 1D directly to the virtual pressure wave $p_{virt}$ ⁠, and the “measured” pressure signal $p$ ⁠, which was calculated from Eq. (7), where 1% of Gaussian noise was added (SNR = 100). For image reconstruction, the inverse matrix was calculated by using the truncated singular value decomposition (T-SVD).²⁵ The SVD of the matrix $M_{p}$ is a factorization of the form $U Σ V^{t}$ ⁠, where $U$ and $V$ are unitary matrices and $Σ$ is a diagonal matrix with non-negative diagonal elements in a decreasing order, called the singular values. For the pseudo-inverse matrix $Σ^{+}$ ⁠, the inverse of the non-vanishing diagonal elements in $Σ$ is taken and the non-significant singular values were set to zero. In the truncated SVD reconstruction, if the singular values get less than 1/SNR, they are set to zero. The pseudo-inverse of the matrix $M_{p}$ is $M_{p}^{+} = V Σ^{+} U^{t}$ ⁠. Then, the reconstructed virtual pressure wave is

p_{rec} = M_{p}^{+} p,

(25)

which is shown in Fig. 2 around the depth of 1 mm and around the depth of 10 mm.

FIG. 2.

View large Download slide

Two layers at a depth $z$ of 1 mm (a) and 10 mm (b), with a layer thickness of 0.2 mm. At a depth of 1 mm, the attenuation of the measured pressure is significantly lower. Therefore, also the reconstructed pressure $p_{rec}$ shows a better spatial resolution. The theoretical resolution given in Table I is $22 μ m$ around 1 mm depth and $71 μ m$ around the depth of 10 mm.

TABLE I.

Truncation frequency and spatial resolution for glycerine.

Propagation distance	Truncation frequency	Spatial resolution
z (mm)	ω_cut/(2π) (MHz)	δ_resolution (μm)
1	42.9	22.4
5	19.2	50.1
10	13.6	70.9

Propagation distance	Truncation frequency	Spatial resolution
z (mm)	ω_cut/(2π) (MHz)	δ_resolution (μm)
1	42.9	22.4
5	19.2	50.1
10	13.6	70.9

B. Virtual wave concept for thermographic reconstruction

In the frequency domain, we can calculate the decay in amplitude of a thermal wave with frequency $ω$ after a distance $z$ either by directly solving the Helmholtz equation (9) or by inserting the Ansatz⁴⁴

T (z, t) = ℜ (T_{0} \exp (i (σ (ω) z - ω t))),

(26)

with the complex wave number $σ (ω) = \sqrt{\frac{i ω}{α}} \equiv \frac{1 + i}{μ}$ and a thin absorbing layer $T_{0} (z) = T_{0} δ (z)$ into the heat diffusion equation (5), which results in

T (z, t) = ℜ (T_{0} \exp (- \frac{z}{μ}) \exp (i \frac{z}{μ} - i ω t)),

(27)

where $ℜ$ is the real part and $μ (ω) \equiv \sqrt{2 α / ω}$ is defined as the thermal diffusion length.⁴⁵ The amplitude of the thermal wave is reduced by a factor $1 / e$ after propagation of that length. The wavenumber or spatial frequency is $k (ω) \equiv 1 / μ (ω) = \sqrt{ω / 2 α}$ ⁠. Similar as for acoustic attenuation described in Sec. II A, for frequencies larger than the truncation frequency $ω_{cut}$ ⁠, the amplitude of these wave components is damped below the noise level, with the truncation wavenumber $k_{cut}$ ⁠, and for the truncation diffusion length, frequency, or wavenumber, we get

\begin{aligned} SNR & \exp (- \frac{z}{μ_{cut}}) = 1, or \\ ω_{cut} & = 2 α {(\frac{\ln (SNR)}{z})}^{2}, or \\ k_{cut} & = \frac{\ln (S N R)}{z} . \end{aligned}

(28)

Similar to Eq. (23), the spatial resolution

δ_{resolution} = \frac{π}{k_{cut}} = \frac{π z}{\ln (SNR)}

(29)

is half of the wavelength of the truncation frequency. For thermographic reconstruction, the resolution limit decays linear with depth $z$ ⁠. In comparison, for photoacoustic reconstruction, the resolution limit decays proportional to the square root of $z$ [Eq. (23)].

The objective of the solution of the inverse problem in thermographic reconstruction is to calculate the virtual wave temperature from the noisy temperature measurements. The one-dimensional solution of the wave equation is a wave package of constant shape traveling with sound velocity, because the group velocity is equal to the phase velocity. Therefore, in 1D the virtual wave immediately gives the reconstruction at time $t = 0$ by projecting it back at a distance $z = c t$ ⁠. Similar to photoacoustic imaging, the optimum reconstruction can be achieved in the temporal frequency domain (⁠ $ω$ -space) by considering frequencies up to the truncation frequency $ω_{cut}$ ⁠. If $T_{virt}$ has the shape of a sharp pulse, ideally a Dirac delta distribution, the pulse will keep its shape over a travel distance, as the virtual wave equation (8) has only solutions without dispersion and attenuation. In $ω$ -space, for ${\tilde{T}}_{virt}$ also the high frequencies will contribute to the signal after some propagation. From the delta pulse in $T_{virt}$ ⁠, the related temperature signal $T$ gets broadened with increasing distance $z$ according to Eq. (19).

In Fig. 3(a), the time-dependent temperature signal calculated using Eq. (21) is shown for three different heat sources at layer depths of 1, 3, and 5 mm. The layer thickness is the same as in the photoacoustic example with 0.2 mm. A low carbon steel (⁠ $< 0.4 % C$ ⁠) with the following material parameters is theoretically investigated: $k = 43 W / (mK)$ ⁠, $C_{p} = 465 J / (kgK)$ ⁠, and $ρ = 7850 {kg / m}^{3}$ ⁠. The resulting thermal diffusivity, which describes the speed of heat propagation, is $α = k / (ρ C_{p}) = 11.8 10^{- 6} m^{2} / s$ ⁠. The characteristic time (or timescale) to reach the maximum temperature in an 1D heat diffusive process is given by $t_{D} = z^{2} / (2 α)$ ⁠. Due to this quadratic dependence, the temperature curves in Fig. 3(a) show very widespread timescales for different layer depths $z$ ⁠. This can be interpreted by the concept of “thermal waves”: the short time behavior is dominated by the fastest propagating, high frequency components of the pulse and the long time behavior by the low frequency components, with a phase velocity of $\sqrt{2 α ω}$ ⁠.⁴⁴ The sample can be considered a low pass filter, because the high frequency components are heavily attenuated as they propagate, and so only the slow moving low frequency components reach the surface from deeper layer sources in the sample. In Fig. 3(b), the reconstructed signal $T_{rec}$ is shown for the three rectangular signals of $T_{virt}$ ⁠, whereby for the regularization the truncated-SVD is used with a truncation value for the smallest singular value of 1/SNR with a SNR of 100. The truncation frequencies are listed in Table II for the three different layer depths. As the resolution limit is the size of a small, or ideally point-like, reconstructed source at a certain depth, the resolution after some propagation distance $z = c t^{'}$ is half the wavelength at this frequency or the width of the reconstructed signal calculated by the T-SVD. This is in good agreement with the spatial resolution in Eq. (29). The thermographic resolution limit is proportional to the travel distance $z$ of the signal between the excitation location and the detection surface and inversely proportional to the natural logarithm of the SNR, whereby the thermal diffusivity has no influence. The resolution limits for a SNR of 100 are given in Table II.

FIG. 3.

View large Download slide

Thermal transients resulting from pulse heating of three different layer sources at a depth $z$ of 1 mm, 3 mm, and 5 mm, with a layer thickness of 0.2 mm and the thermal properties of stainless steel. The reconstructed temperature signals from the virtual signals in three different depths and the decreasing resolution with traveling distance $z = c t^{'}$ are shown in (b). The resolution is in the order of the layer depth.

TABLE II.

Truncation frequency and spatial resolution for steel in thermographic reconstruction with a SNR = 100.

Propagation distance	Truncation frequency	Spatial resolution
z (mm)	ω_cut/(2π) (Hz)	δ_resolution (mm)
1	79.5	0.68
3	8.8	2.05
5	3.2	3.41

Propagation distance	Truncation frequency	Spatial resolution
z (mm)	ω_cut/(2π) (Hz)	δ_resolution (mm)
1	79.5	0.68
3	8.8	2.05
5	3.2	3.41

III. THERMODYNAMICS AND INFORMATION

The process of information gaining during photoacoustic and photothermal imaging is a physical one which has to obey the laws of (non-equilibrium) thermodynamics. Determining the virtual pressure or temperature wave from the measured pressure or temperature is an ill-posed or an ill-conditioned inverse problem, which needs regularization.⁴ In Sec. II, the truncation frequency $ω_{cut}$ was chosen as a regularization parameter, where it was assumed by regularization using the T-SVD method that all signal components below $ω_{cut}$ could be used for reconstruction, and all components above $ω_{cut}$ ⁠, where the signal is damped below the noise level, provide no additional information for reconstruction. Therefore, the signal amplitude compared to the noise level plays an important role in determining the truncation frequency $ω_{cut}$ ⁠, which serves as a regularization parameter for the inverse problem. Where does this noise come from and what is its relation to information gaining in imaging?

Before answering these essential questions, some basics of thermodynamics and statistical physics should be reviewed. A microstate is a specific microscopic configuration of a thermodynamic system, e.g., with certain locations or velocities of all of the molecules. A microstate can be represented as a single point with coordinates $x$ in a usually very high dimensional phase space. In contrast, the macrostate of a system refers to a few macroscopic properties, such as its temperature, pressure, volume, or density. A macrostate is characterized by a probability distribution $ρ (x)$ in the phase space of possible microstates. This distribution describes the probability of finding the system in a certain microstate. In classical mechanics, the position and momentum variables of a particle can vary continuously, so the set of microstates is actually uncountable and $x$ gets a continuous variable. The time evolution of all the particles is determined by the Hamiltonian $H$ ⁠, which gives the total energy of a microstate, the sum of the kinetic and potential energy. The systems energy

U_{t} = \int H (x) ρ_{t} (x) d x

(30)

at time $t$ is the mean value or also called expectation value of the system Hamiltonian $H$ ⁠, where $ρ_{t} (x)$ is the probability distribution at time $t$ ⁠.

Here, we follow the work of Esposito and van den Broeck⁷ about the second law and the connection between entropy production and information loss. They have used a general Hamiltonian system in contact with an ideal heat bath at temperature $T$ ⁠. $W_{t}$ is the work performed on the system, and $Q_{t}$ is the heat coming from the system into the ideal heat bath. In our case, the system is the sample, the performed work is the energy $W$ of the short laser pulse at the time $t = 0$ ⁠, and $Q_{t}$ is the diffused heat or the dissipated heat from the acoustic attenuation until time $t$ ⁠. The ideal heat bath is the samples environment having an ambient temperature $T$ (see Fig. 4). The system entropy is the Shannon entropy

S_{t} = - k_{B} \int ρ_{t} (x) \ln (ρ_{t} (x)) d x,

(31)

where $k_{B}$ is the Boltzmann constant, and the corresponding system free energy is given by

F_{t} = U_{t} - T S_{t} .

(32)

Here, $T$ is the temperature in Kelvin of an ideal heat bath with which the system is in contact—the ambient temperature of the sample. The excitation laser pulse deposits the work $W$ at the time $t = 0$ ⁠, which either diffuses as heat out of the system or is the heat from dissipation of the photoacoustic pulse by acoustic attenuation, called $Q_{t}$ from the laser pulse until the time $t$ ⁠. After a long time, the system is in the same thermal equilibrium as in the time before the excitation pulse, and all the energy $W$ of the laser pulse has left the system and is in the heat bath (⁠ $Q_{t} = W$ ⁠), as sketched in Fig. 4.

FIG. 4.

View large Download slide

Illustration of the thermodynamics of the imaging process in phase space: a system in equilibrium state $ρ_{eq}$ is kicked at a time $t = 0$ with magnitude $x_{0}$ and performed work $W$ to a state $ρ_{kick}$ far from equilibrium, followed by a dissipative or diffusive process back to equilibrium. Until the time $t$ ⁠, the heat $Q_{t}$ flows to the surroundings of the system at an ambient temperature $T$ ⁠. At $t \to \infty$ ⁠, equilibrium is reached again and $Q_{t}$ gets $W$ ⁠. $x$ is the coordinate in phase space or a set of reduced variables, which captures the information on the work (see the text). The arrows connecting $ρ_{kick}$ at time $t = 0$ ⁠, $ρ_{t}$ at $t > 0$ ⁠, and $ρ_{eq}$ at $t \to \infty$ indicate the tube of trajectories, which is “thin” for macroscopic systems as deviations from the mean values $x (t)$ are small.

Following conservation of energy (first law of thermodynamics), the corresponding energy change $Δ U_{t} := U_{t} - U_{eq}$ of the system is given by

Δ U_{t} = W - Q_{t} for t > 0.

(33)

At $t \to \infty$ ⁠, equilibrium is reached again, $Q_{t}$ gets $W$ ⁠, and the change in system energy $Δ U_{t} = 0$ ⁠.

If the system is in equilibrium with the surroundings at the ambient temperature $T$ ⁠, the equilibrium distribution is the canonical distribution

\begin{aligned} ρ_{eq} (x) = \frac{1}{Z} \exp (- β H (x)) \\ with \\ Z = \int \exp (- β H (x)) d x, \end{aligned}

(34)

where the normalization function of the equilibrium distribution is called the canonical participation function $Z$ and the thermodynamic $β = 1 / (k_{B} T)$ ⁠. This Maxwell-Boltzmann distribution $ρ_{eq}$ can be determined by minimizing the entropy $S_{t}$ under the two constraints on the total particle number and on the average energy per particle, e.g., Refs. 46 and 47. For the equilibrium entropy, one gets

S_{eq} = k_{B} \int ρ_{eq} (x) (\ln Z + β H (x)) d x = k_{B} \ln Z + \frac{U_{eq}}{T},

(35)

using Eqs. (31) and (30), and for the equilibrium free energy

F_{eq} = U_{eq} - T S_{eq} = - k_{B} T \ln Z .

(36)

According to the second law of thermodynamics, the entropy of an adiabatically insulated system increases monotonically until thermodynamic equilibrium is established, where the relative entropy $- k_{B} D (ρ_{t} ∣ ρ_{eq})$ gets zero, with the Kullback–Leibler divergence

D (ρ_{t} | ρ_{eq}) = \int ρ_{t} (x) \ln (\frac{ρ_{t} (x)}{ρ_{eq} (x)}) d x

(37)

between the non-equilibrium distribution $ρ_{t}$ at time $t$ and the equilibrium distribution $ρ_{eq}$ ⁠.³⁴ The Chernoff–Stein’s Lemma gives a precise meaning to $D (f ∣ g)$ as a “distance” between two distributions: if $n$ data from $g$ are given, the probability of guessing incorrectly that the data come from $f$ is bound by the error $ϵ = \exp (- n D (f ∣ g))$ ⁠, for $n$ large.^34,48 The free energy of a non-equilibrium state $ρ_{t}$ is higher than that of the equilibrium state by an amount equal to the temperature times the information $I_{t}$ needed to specify the non-equilibrium state

\begin{aligned} T I_{t} \equiv k_{B} T D (ρ_{t} | ρ_{eq}) & = - T S_{t} + \int H (x) ρ_{t} (x) d x + k_{B} T \ln Z \\ = U_{t} - T S_{t} + k_{B} T \ln Z \\ = F_{t} - F_{eq} \equiv Δ F_{t} . \end{aligned}

(38)

For the first equality, Eqs. (37), (31), and (34) are used, and the second equation uses Eq. (30). The third equation uses the definition of the free energy in Eq. (32) and for the equilibrium state in Eq. (36). Due to the excitation pulse at time $t = 0$ ⁠, the free energy jumps from $F_{eq}$ to $F_{kick} = F_{eq} + W$ ⁠, because the pulse shifts the distribution to a higher energy $Δ U_{kick} = W$ [Eq. (33)] but does not change the entropy (Fig. 4), $S_{kick} = S_{eq}$ ⁠.

The change in the non-equilibrium system entropy $Δ S_{t} = S_{t} - S_{eq}$ can be written as a reversible contribution to the heat flow $- Q_{t} / T$ and of an irreversible contribution called entropy production

Δ S_{t} = Δ_{i} S_{t} - \frac{Q_{t}}{T} .

(39)

Together with the definition of the free energy in Eq. (32), the first law equations (33) and (38) we obtain for the entropy production

\begin{aligned} T Δ_{i} S_{t} & = T Δ S_{t} + Q_{t} = Δ U_{t} - Δ F_{t} + Q_{t} \\ = W - Δ F_{t} = W - T I_{t} for t > 0. \end{aligned}

(40)

The free energy $Δ F_{t}$ jumps at time $t = 0$ from zero to $W$ and then decreases and gets zero, the information according to Eq. (38) makes the jump to $W / T$ and then decreases to zero. $Δ I_{t}$ is defined as the information loss, and according to Eq. (40), we get the very important result that the entropy production,

Δ_{i} S_{t} = Δ I_{t} = W / T - I_{t},

(41)

is equal to the information loss.

For a macroscopic sample, the fluctuations that are the variance or the “noise” around the mean values $x (t)$ are small and the change of the “shape” of the distribution during time evolution shown in Fig. 4 can be neglected. Therefore, the change of the Shannon system entropy can be neglected [ $Δ S_{t} \approx 0$ in Eq. (39)] and the change in free energy is equal to the change in system energy, $Δ F_{t} \approx Δ U_{t}$ ⁠, which gives from Eq. (38) for the information

\begin{aligned} I_{t} & = k_{B} D (ρ_{t} ∣ ρ_{eq}) = \frac{Δ F_{t}}{T} \approx \frac{Δ U_{t}}{T} \\ = \frac{W - Q_{t}}{T} \approx k_{B} D (ρ_{kick} ∣ ρ_{eq}) - k_{B} D (ρ_{kick} ∣ ρ_{t}) . \end{aligned}

(42)

In this very good approximation for macroscopic samples, the Kullback–Leibler divergence gets a distance in the mathematical sense, and the energy $W - Q_{t}$ ⁠, which has not yet dissipated or diffused to the surroundings at temperature $T$ ⁠, divided by this temperature is the information content. After a long time, when all the work $W$ has dissipated or diffused, $I_{t}$ gets zero and the equilibrium distribution $ρ_{eq}$ is reached. But according to Chernoff–Stein’s Lemma after some time $t_{cut}$ ⁠, the distribution $ρ_{t_{cut}}$ cannot be statistically distinguished from $ρ_{eq}$ ⁠, if $k_{B} T D (ρ_{t} ∣ ρ_{eq})$ gets smaller than the equilibrium energy $U_{eq}$ ⁠.

For real world examples having a high dimensional phase space $x$ ⁠, the actual distribution density $ρ_{t}$ can be hardly determined, but for the information loss from Eq. (42) it is sufficient to evaluate the dissipated or diffused work from the mean-value equations. Usually, a set of reduced variables with a drastically lower dimensionality than the phase space can be used, which captures the information on the mean dissipated or diffused work. In the macroscopic approximation, the actual noise or fluctuation around the mean value is irrelevant. Only its amplitude is significant for setting the lower bound in Chernoff–Stein’s Lemma. Subsections III A and III B give a small example for dissipation and for diffusion to clarify the relation between entropy production and information loss.

A. Kicked Brownian particle

We now consider the velocity $v$ of a particle as a stochastic process, which was used to describe the Brownian motion of a particle, but for simplicity only one velocity component (1D) is considered here. The essential relation between fluctuation (noise) and dissipation (viscous damping) is shown. Stochastic processes can be described mathematically, e.g., by Master equations, Langevin or Fokker–Planck equations shown, e.g., in the books of van Kampen,⁴⁹ Gardiner,⁵⁰ or Risken.⁵¹ In the Langevin equation

\frac{d v (t)}{d t} = - γ v (t) + σ η (t),

(43)

the environmental forces on a particle in Newton’s law are a linear damping term together with random noise. The linear damping $- γ v$ is a viscous drag, and $σ$ is the amplitude of the white noise $η$ ⁠, which has a zero mean value and is uncorrelated in time: $⟨ η (t) η (t^{'}) ⟩ = δ (t - t^{'})$ ⁠. The Langevin equation governs an Ornstein–Uhlenbeck (O-U) process, named after L. S. Ornstein and G. E. Uhlenbeck, who formalized the properties of this continuous Markov process.⁵² It was shown by using the statistical properties and the continuum limit of the white noise $η$ that the Langevin equation is equivalent to a description based on a Fokker–Planck equation. For the time-dependent distribution density $ρ_{t} (v)$ of the velocity, a linear differential equation

\frac{\partial ρ_{t} (v)}{\partial t} = \frac{\partial (γ v ρ_{t} (v))}{\partial v} + \frac{σ^{2}}{2} \frac{\partial^{2} ρ_{t} (v)}{\partial v^{2}}

(44)

can be derived.^51,53

We start with the equilibrium state [zero time derivative if inserted into Eq. (44)] as the initial velocity distribution

ρ_{eq} (v) = \frac{1}{Z} \exp (- \frac{γ}{σ^{2}} v^{2}) with Z = σ \sqrt{\frac{π}{γ}},

(45)

which is a Gaussian distribution with mean value zero and a variance $Var (v) = σ^{2} / (2 γ)$ ⁠, representing the statistical “noise.” Comparison with Eq. (34) and the kinetic energy $H (v) = m v^{2} / 2$ gives

β H (v) = \frac{1}{k_{B} T} \frac{m v^{2}}{2} = \frac{γ}{σ^{2}} v^{2} or Var (v) = \frac{σ^{2}}{2 γ} = \frac{k_{B} T}{m} .

(46)

This relation states a connection between the strength of the fluctuations, given by $σ$ ⁠, and the strength of the dissipation $γ$ ⁠. This is the fluctuation–dissipation theorem for uncorrelated white noise.

At a time zero, the particle is kicked, which causes an immediate change in velocity of $v_{0}$ (kick magnitude) and the distribution density after the kick is

ρ_{kick} (v) = ρ_{eq} (v - v_{0}) = \frac{1}{Z} \exp (- \frac{γ}{σ^{2}} {(v - v_{0})}^{2}) .

(47)

The solution of Eq. (44) for the time-dependent distribution density $ρ_{t} (v)$ at $t > 0$ with $ρ_{kick} (v)$ as initial condition at $t = 0$ is

\begin{aligned} ρ_{t} (v) & = \frac{1}{Z} \exp (- \frac{γ}{σ^{2}} {(v - \bar{v} (t))}^{2}) \\ with \\ \bar{v} (t) & = v_{0} \exp (- γ t) for t > 0, \end{aligned}

(48)

which gives a Gaussian distribution with time dependent mean value $\bar{v} (t)$ but a constant variance $Var (v) = σ^{2} / (2 γ)$ ⁠. As shown by Burgholzer and Hendorfer, this is a general feature of Gauss–Markov processes, also for higher dimensions: taking the equilibrium as an initial condition results for all times after the kick in a distribution with a constant (co)variance (matrix) equal to the equilibrium variance.⁹ One realization of the kicked O-U process is shown in Fig. 5. At the time $t = 0$ ⁠, a kick with a magnitude of $v_{0} = 10$ occurs. For a time $t > 0$ ⁠, the information about the magnitude of the kick gets more and more lost due to the fluctuations. Now, this increasing information loss is quantified and compared to the mean entropy production.

FIG. 5.

View large Download slide

The circles show a typical realization of a kicked Ornstein–Uhlenbeck process defined by the Langevin equation (43). The scaled time $t γ$ is on the horizontal axis. The velocity $v$ on the vertical axis is scaled to have a unit variance. At the time $t = 0$ ⁠, a kick magnitude of $v_{0} = 10$ is added to the scaled velocity. Increasing in time the information about the magnitude of the kick gets more and more lost due to the fluctuations. The solid line represents the mean, and the dashed lines represent the mean $\pm$ standard deviation, which is the square root of the variance.

According to Eqs. (30) and (33) due to the kick at time zero, the energy of the particle is increased by $W = m v_{0}^{2} / 2$ and

\begin{aligned} Δ U_{t} & = \int H (v) ρ_{t} (v) d v - \int H (v) ρ_{eq} (v) d v \\ = \frac{m \bar{v} {(t)}^{2}}{2} = W \exp (- 2 γ t) . \end{aligned}

(49)

From Eq. (31), the entropy

S_{t} = - k_{B} \int ρ_{t} (v) \ln (ρ_{t} (v)) d v = k_{B} \ln Z + \frac{k_{B}}{2}

(50)

stays constant in time (⁠ $Δ S_{t} = 0$ ⁠), as the distribution is shifted in time but with a constant “shape.” Therefore, the approximation in Eq. (42) gets exact and the information content is

\begin{aligned} I_{t} & = k_{B} D (ρ_{t} | ρ_{eq}) = k_{B} \int ρ_{t} (v) \ln (ρ_{t} (v) / ρ_{eq} (v)) d v \\ = k_{B} \frac{γ}{σ^{2}} \bar{v} {(t)}^{2} = \frac{k_{B}}{2} \frac{v_{0}^{2}}{Var (v)} \exp (- 2 γ t) . \end{aligned}

(51)

According to Chernoff–Stein’s Lemma, if the information $I_{t}$ at a time $t_{cut}$ gets less than a certain limit, here chosen as the zero-point energy divided by the temperature (⁠ $U_{eq} / T = k_{B} / 2$ ⁠), the distribution $ρ_{t} (v)$ cannot be statistically distinguished from the equilibrium distribution $ρ_{eq} (v)$ ⁠. The information content is

\begin{aligned} I_{t_{cut}} = k_{B} D (ρ_{t_{cut}} | ρ_{eq}) & = \frac{k_{B}}{2} \frac{v_{0}^{2}}{Var (v)} \exp (- 2 γ t_{cut}) \\ = \frac{U_{eq}}{T} = \frac{k_{B}}{2}, \end{aligned}

(52)

which gives

t_{cut} γ = \ln (\frac{v_{0}}{\sqrt{Var (v)}}) \equiv \ln (SNR) .

(53)

At this truncation time $t_{cut}$ ⁠, not only the distribution cannot be distinguished from equilibrium, but also the mean velocity $\bar{v} (t) = v_{0} \exp (- γ t)$ gets less than the noise level $v_{0} / SNR$ ⁠. In Fig. 5, this is at the scaled time $t_{cut} γ = \ln (10) \approx 2.3$ ⁠. This coincidence, that the information related criterion in Eq. (52), gives the same time $t_{cut}$ when the signal gets less than the noise always happens, if the energy described by the Hamiltonian $H$ is proportional to the square of the signal amplitude.

B. Ideal gas diffusion between two boxes

Another example for an Ornstein–Uhlenbeck process using diffusion instead of dissipation is the distribution of $2 N$ particles of an ideal gas between two boxes of volume $V$ ⁠, which are connected by a small hole where particles can slip through between the two boxes (effusion). The constant time rate $γ$ is the reciprocal value of the mean time of a particle to stay in one of the boxes before it slips through the hole. At a certain time $t$ ⁠, $N + δ N (t)$ particles are situated in the left box and $N - δ N (t)$ particles are in the right box (Fig. 6).

FIG. 6.

View large Download slide

Distribution of $2 N$ particles between two boxes of volume $V$ ⁠, which are connected by a small hole where particles can change place between the two boxes at a constant time rate $γ$ (effusion). At a certain time $t$ ⁠, $N + δ N (t)$ particles are situated in the left box and $N - δ N (t)$ particles are in the right box.

$δ N (t)$ can be approximated by an Ornstein–Uhlenbeck process for a big number of particles (see, e.g., the tutorial introduction to stochastic processes by Lemons and Langevin⁵⁴). At the time $t = 0$ ⁠, $N_{0}$ particles are “pumped” from the right to the left box. The mean value $δ \bar{N} (t) = N_{0} \exp (- γ t)$ for $t > 0$ is given by an exponential decay in time, and the variance of the Gauss distribution is again constant in time with $Var (δ N) = \frac{σ^{2}}{2 γ}$ ⁠, according to Eq. (48). However, $Var (δ N)$ is not an energy, and therefore $σ^{2}$ cannot be expressed in terms of a rate $γ$ and a temperature T by requiring the equipartition of energy at equilibrium. The fluctuation–dissipation theorem as in Eq. (46) does not apply; $Var (δ N)$ fluctuates without dissipation.⁵⁴ The information loss can be described analogously to Eq. (51) by an exponential decay

I_{t} = k_{B} \frac{γ}{σ^{2}} δ \bar{N} {(t)}^{2} = \frac{k_{B}}{2} \frac{N_{0}^{2}}{Var (δ N)} \exp (- 2 γ t) .

(54)

Spatial diffusion is the cause for the entropy production and is described according to Boltzmann by $S = k_{B} \ln W$ ⁠. $W$ is the number of possibilities to distribute the particles into the two boxes for a certain $δ \bar{N}$ ⁠. Using Stirlings formula for a big number of particles $N$ ⁠, one gets for the entropy⁴⁶

Δ S_{t} = - k_{B} \frac{1}{N} δ \bar{N} {(t)}^{2} = - I_{t} .

(55)

By comparison with Eq. (54), one gets again a relation between the fluctuations, given by $σ$ ⁠, and the strength of the diffusion, given by the rate $γ$ (fluctuation–“diffusion” theorem) and for the variance of $δ N$ ⁠, we get

Var (δ N) = \frac{σ^{2}}{2 γ} = \frac{N}{2} .

(56)

IV. IMAGING AS AN INVERSE PROBLEM

In Sec. III, imaging with an excitation pulse at time $t = 0$ has been described by means of nonequilibrium statistical physics. Pressure or temperature is proportional to the momentum or the kinetic energy of many particles, which fluctuate around their mean value. For macroscopic samples, these fluctuations are very small and can usually be neglected in the thermodynamic limit. But for inverse problems, these fluctuations are highly “amplified,” which was shown in Sec. III for a system kicked out of the equilibrium, followed by a dissipative process back to equilibrium (Fig. 4). The inverse problem of estimating the kick-magnitude from an intermediate state a certain time after the kick is ill-posed. Just after the kick, its magnitude can be estimated very well. A long time after the kick the state is nearly in equilibrium and all the information about the kick magnitude is lost. For macroscopic systems, it could be shown that the information loss is in a good approximation just the mean dissipated energy or diffused heat divided by the temperature, which is the mean entropy production. In imaging, the spatial resolution and the information content are strongly correlated, and therefore a loss of information results in a loss of resolution, which is quantified for 1D examples in Secs. II A and II B.

In Secs. III A and III B, the mean dissipated energy and the resulting loss of information are calculated explicitly for a kicked Brownian particle and for an ideal gas diffusing between two boxes. For these study cases, also the time-dependent probability distribution could be calculated explicitly, and the deduced information losses are in agreement with the results from mean-value calculations, even when the distributions are broad and still far away from the thermodynamic limit. A truncation time could be given, for which the state after that time cannot be distinguished from the equilibrium distribution according to Chernoff–Stein’s Lemma, when the Kullback–Leibler divergence gets too small. The amplitude of the fluctuations and the dissipation or the diffusion are not independent. Fluctuations and mean entropy production can be thought as two sides of the same coin and are connected by the fluctuation–dissipation theorem. For systems near thermal equilibrium in the linear regime, such relations between entropy production and fluctuation properties have been found by Callen and Greene,⁵⁵ Callen and Welton,⁵⁶ and Greene and Callen.⁵⁷ This fluctuation–dissipation theorem is a generalization of the famous Johnson⁵⁸ and Nyquist⁵⁹ formula in the theory of electric noise. It is based on the fact that in the linear regime the fluctuations decay according to the same law as a deviation from equilibrium following an external perturbation.

For macroscopic systems, it is not necessary to describe the full stochastic process to get the influence of the fluctuations. The relevant information loss can be calculated from the averaged behavior (mean value), which describes the usually known macroscopic evolution of the system in time, as the mean dissipated work or diffused heat divided by the temperature. This remarkable feature might be the reason that regularization methods for ill-posed inverse problems work so well, although they use only the mean-value equations and not the detailed stochastic process to describe the time evolution. The choice of an adequate regularization parameter, e.g., the truncation value for the truncated singular value decomposition (SVD) method, is equivalent to the choice of the error level $ϵ$ in Chernoff–Stein’s Lemma.

A prominent class of ill-posed inverse problems is non-destructive imaging (NDI), where the information about the spatial pattern of a sample’s interior has to be transferred to the sample surface by certain waves, e.g., ultrasound or thermal waves (Fig. 1). Imaging is done by reconstruction of the interior structure from the signals measured on the sample surface, e.g., by backprojection or time reversal for a photo-acoustically induced ultrasound wave.^11,60 There are several effects that limit the spatial resolution for photoacoustic imaging. Beside insufficient instrumentation and data processing, one principle limitation comes from attenuation of the acoustic wave or from heat diffusion (Sec. II). Attenuation during wave propagation can be caused by dissipation of acoustic energy to heat or by acoustic scattering. In this article, it is assumed that, as with dissipation, the information is also lost with scattering and that scattered waves are not used for image reconstruction. Using scattered waves for reconstruction would be possible, but this is beyond the scope of this Tutorial. Therefore, the entire information loss due to acoustic attenuation is lost and cannot be compensated by subsequent processing algorithms. In this section, it is shown that the loss of information, which is equal to the entropy production (Sec. III) as the dissipated energy or diffused heat divided by the temperature, is a principle thermodynamic limit, which cannot be compensated. Using the information loss and entropy production for a kicked process derived in Sec. III, it is shown that the resolution limit depends just on the macroscopic mean-value equations and is independent of the actual stochastic process, as long as the macroscopic equations describe the mean work and therefore also the mean dissipated work.

In this section, we do not attempt to model the process of acoustic attenuation or heat diffusion as a stochastic process, which we have done earlier as a Gaussian process.⁹ We use from Sec. II the mean-value equations for the pressure $p (r, t)$ or temperature $T (r, t)$ evaluation in the frequency domain, $\tilde{p} (r, ω)$ or $\tilde{T} (r, ω)$ ⁠, respectively.

For an acoustic wave in the frequency domain, according to Eq. (10), the amplitude of the wave component with frequency $ω$ is damped by the factor $\exp (- α_{0} ω^{2} r)$ after propagating a distance $r$ ⁠. The energy $Δ U_{ω}$ of the acoustic wave with frequency $ω$ is proportional to the square of the pressure amplitude

Δ U_{ω} = 0.5 χ Δ V {| \tilde{p} (r, ω) |}^{2} = 0.5 χ Δ V \exp (- 2 α_{0} ω^{2} r),

(57)

which can be found, e.g., in Morse and Ingard,⁶¹ where $χ = 1 / (c^{2} ρ)$ is the adiabatic compressibility with the density $ρ$ and $Δ V$ is the measurement volume. Using that $Var (ρ) = k_{B} T / (χ Δ V)$ is the variance of the pressure (e.g., from Landau and Lifshitz⁶²) one gets from Eq. (42) for the information content of the Fourier component with frequency $ω$ ⁠: $I_{ω} = Δ U_{ω} / T = 0.5 k_{B} {SNR}^{2} \exp (- 2 α_{0} ω^{2} r)$ ⁠. The signal-to-noise ratio at the distance $r = 0$ is the reciprocal value of the square root of the variance of the pressure, as the signal amplitude in the frequency domain is normalized to one [Eq. (10)]. The truncation frequency $ω_{cut}$ is defined in the manner that the information content for frequencies larger than $ω_{cut}$ is so low that its distribution cannot be distinguished from the equilibrium distribution within a certain statistical error level (Chernoff–Stein’s lemma). This level is $Δ U_{eq} / T = 0.5 χ Δ V Var (p) / T = 0.5 χ Δ V k_{B} T / (χ Δ V) / T = 0.5 k_{B}$ ⁠, which results in

I_{ω_{cut}} = 0.5 k_{B} {SNR}^{2} \exp (- 2 α_{0} ω_{cut}^{2} r) = 0.5 k_{B} .

(58)

This gives the same truncation frequency as in Eq. (22), where the amplitude is damped just below the noise level. As described already in Sec. III A, the information related criterion in Eq. (58) gives the same truncation value as when the signal amplitude gets less than the noise, if the energy $Δ U_{ω}$ is proportional to the square of the signal amplitude $| \tilde{p} (r, ω) |$ [Eq. (57)].

Similarly, also for particle diffusion, the information content $I_{t}$ is proportional to the square of the mean difference in the particle number $δ \bar{N} (t)$ [Eq. (54)]. For heat diffusion and thermal waves, the information content is proportional to the square of the temperature deviation from the equilibrium temperature, in real time space and in the frequency domain. Heat diffusion can be described as an Ornstein–Uhlenbeck process.^9,63 Therefore, the information related criterion in Eq. (42) for thermal waves gives the same truncation frequency $ω_{cut}$ when the signal gets less than the noise with the results given in Sec. II B.

As mentioned at the end of Sec. II, imaging in two (2D) or three dimensions (3D) can always be reduced using a two-stage process: first, for each detector location the virtual pressure signal in the absence of attenuation or the virtual thermal wave is calculated from the measured acoustic or thermal signal, respectively. This is a one-dimensional (1D) reconstruction problem, described in discrete time steps by Eqs. (17) and (21). In a second step, any reconstruction method for photoacoustic tomography without acoustic attenuation, such as time reversal or backprojection, can be used for reconstructions in higher dimensions.¹¹ In 1D, the virtual wave immediately gives the reconstruction at time $t = 0$ by projecting it back at a distance $c t$ (see Secs. II A and II B). Therefore, it is sufficient to examine the acoustic attenuation or thermal diffusion of 1D waves and the reconstruction in 1D. Compensation of acoustic attenuation or thermal diffusion in higher dimensions can always be reduced to 1D, which was also explicitly shown for signals from a layer (1D), cylinder (2D), and a sphere (3D).²⁰ Therefore, in the beginning, we will show a 1D example.

A. 1D imaging

1. Acoustic attenuation in a fat tissue layer

For excitation of plane acoustic waves, the surface of a silicon wafer was illuminated by a nanosecond laser pulse (Fig. 7). Due to the abrupt local heating due to absorption of the 532 nm wavelength laser light pulse, the silicon wafer expands thermoelastically and generates an acoustic wave, propagating through a layer of fatty tissue and is detected by an unfocused ultrasound transducer (V358-SU, Panametrics, Waltham, MA). The measured acoustic pressure without fat, through 6 mm thick porcine fat tissue and through 20 mm thick tissue, is shown in Fig. 8(a) as a function of time.^64,65 For comparison, the measured signals are scaled and time-shifted. The frequency spectrum was calculated by Fourier transformation in time and the frequency dependent attenuation for the 6 mm and 20 mm fat layer was determined by dividing their spectrum through the spectrum of the reference signal without fat. It turned out that a power law $α (ω) = α_{0} ω^{n}$ with an exponent $n = 1.5$ and $α_{0} = 0.87 dB {MHz}^{- n} {cm}^{- 1}$ fits the attenuation in a wide frequency range very well.⁶⁴

FIG. 7.

View large Download slide

Setup for the generation and detection of acoustic plane waves. Abrupt local heating of a silicon wafer by nanosecond laser pulses leads to the emission of strong broadband ultrasonic plane waves. Porcine subcutaneous fat tissue in the propagation path induces frequency dependent attenuation of the acoustic signals. The fat tissue is fastened between two aperture disks applying a small axial force on the tissue. Distance bolts with 6 mm or 20 mm ensure two precise lengths of the attenuation path. For these two lengths, attenuated acoustic plane waves were detected by an unfocused piezoelectric transducer, which was aligned by worm screws to ensure a one-dimensional signal propagation and detection. Image from Burgholzer et al., J. Imaging 5, 13 (2019). Copyright 2019 Authors(s), licensed under CC BY 4.0.⁶⁶

FIG. 8.

View large Download slide

Acoustic attenuation in a fat tissue layer: (a) Measured acoustic pressure as a function of time without fat tissue (“no fat”), with 6 mm thick porcine fat tissue and with 20 mm porcine fat tissue. To make comparison easier, the signals were time-shifted and scaled and (b) the reconstruction results using T-SVD for regularization to compensate attenuation in fatty tissue of 6 mm thickness and 20 mm thickness. This corresponds to a spatial resolution limit of $32 μ m$ for 6 mm fat and $70 μ m$ for 20 mm fat, resulting from entropy production. The matrix $M_{z}$ was multiplied by the convolution matrix $M_{water}$ of the water signal to get a $δ$ -like pulse for the pure water signal without fatty tissue in the measurement chamber.

The experimental setup now is slightly different compared to the photoacoustic imaging setup described in Sec. II. In photoacoustic imaging, the short excitation pulse at the time zero excites a pressure wave in the whole sample at the same time. If an acoustic signal arrives later at the detector, this comes from a longer propagation distance and causes a higher attenuation. Here, the fat layer always has the same thickness and therefore the attenuation for earlier or later arriving parts of the signal at the detector is the same. Only for the excitation as a Dirac delta in space and time as described in Eq. (16), this is the same. Since the differential equations are linear, the general attenuated signal is a temporal convolution of the input signal with the attenuated solution for Dirac excitation, or a simple multiplication in the frequency domain.¹⁶ The discrete version of the relationship between the measured pressure and the virtual pressure wave now, instead of Eq. (17), becomes

p = M_{z} p_{virt},

(59)

where $p$ and $p_{virt}$ are the vectors of the attenuated and virtual pressure signal at discrete time steps, respectively. Writing the discrete Fourier and inverse Fourier transformation as multiplication by $F$ and its conjugate transpose $F^{*}$ ⁠, and $diag (.)$ forms a diagonal matrix, $M_{z}$ can be written by neglecting the dispersion as⁶⁴

M_{z} = F^{*} diag (\exp (- α_{0} ω^{n} z)) F,

(60)

which immediately shows that the singular values of $M_{z}$ decay exponentially and according to Eq. (22) for larger frequencies than the truncation frequency $ω_{cut}$ the amplitude of these wave components is damped below the noise level, which results in

SNR \exp (- α_{0} ω_{cut}^{n} z) = 1 or ω_{cut} = \sqrt[n]{\frac{\ln (SNR)}{α_{0} z}} .

(61)

The spatial resolution according to Eq. (23) is half the wavelength at the truncation frequency $ω_{cut}$ ⁠. For the SNR of 1358 (63 dB), the truncation frequency for a fat thickness of 6 mm and 20 mm is 24 MHz and 11 MHz, respectively, which corresponds to a spatial resolution of $32 μ m$ after 6 mm of fat and $70 μ m$ after 20 mm of fat.

The acoustic attenuation in water compared to fat can be neglected.⁶⁴ Therefore, the signal measured in water without fat can be taken as the virtual pressure signal $p_{virt}$ ⁠, which is no single Dirac delta pulse, but gets negative and shows additional “ringing” because of the laser ultrasound excitation within the silicon wafer and the characteristic of the piezoelectric transducer and the amplifier. All these influences on the signal can be described by a matrix $M_{water}$ with $p_{virt} = M_{water} p_{δ}$ ⁠, and Eq. (59) can be written as

p = M_{z} M_{water} p_{δ} .

(62)

Using the truncated SVD for the inversion of $M_{z} M_{water}$ ⁠, the reconstructions of $p_{δ}$ from the measurements $p$ for water, 6 mm fat, and 20 mm fat are shown in Fig. 8(b). Multiplying the temporal width of these peaks of 21 ns and 46 ns by the sound velocity of 1512 m/s in fatty tissue perfectly fits to the resolution limits of $32 μ m$ after 6 mm of fat and $70 μ m$ after 20 mm of fat derived above from the truncation frequencies.

2. Axial thermal profiling of planar heat sources

To generate an internal planar heat source for an experimental study, a thin plate (foil) was embedded in pure epoxy resin. Laser light is absorbed for thermal excitation and heats the foil [Fig. 9(b)]. In the case of an electrically conductive film, electromagnetic induction can also be used to heat the film [Fig. 9(a)]. The epoxy resin is opaque in the spectral sensitivity range of the mid-wave IR camera (⁠ $3 - 5 μ m$ ⁠). Therefore, the measured infrared radiation only comes from the sample surface (⁠ $z = 0 mm$ ⁠).

FIG. 9.

View large Download slide

Experimental setup for 1D thermographic reconstruction of internal heat sources: a thin foil (thickness approx. 0.2 mm) embedded in epoxy resin is heated inductively (a) or by the absorption of laser light (b). Since the sample diameter (40 mm) is significantly larger than the deepest defect depth (11.7 mm), we can assume one-dimensional heat conduction in the axis of symmetry.

The thermal diffusivity of the epoxy sample was obtained by the Parker method⁶⁷ as $α = 1.4 \times 10^{- 7} m^{2} s^{- 1}$ ⁠. For two different epoxy samples, with the embedded film at a depth of $z = 7.2 mm$ and 11.7 mm, the surface temperature as a function of time has been measured with a quantum IR detector [Fig. 10(a)]. Due to the low temperature increase on the sample surface, the measurement noise is approximately additive white Gaussian noise with a standard deviation of 25 mK.⁶⁸ The measurement parameters are listed in Table III, with the pulse duration $t_{p}$ ⁠, the number of discretizations in time and space $N$ ⁠, and the temporal and spatial discretization $Δ t$ and $Δ z$ ⁠.

FIG. 10.

View large Download slide

Depth resolved thermal source profiling: (a) Thermal transients resulting from inductive pulse heating for two different layer sources with a layer thickness of 0.2 mm at a depth $z$ of 7.2 mm and 11.7 mm and (b) the solution of regularized inverse problem using the direct method of T-SVD and the iterative procedure ADMM including prior information. The reconstructed depth profiles $T_{rec}$ are compared with the ideal source distribution $T_{virt}$ ⁠.

TABLE III.

Measurement and reconstruction parameters for the 1D reconstruction of the internal heat sources. The SNR was derived from the considered singular values.

z	t_p	N	Δt	Δz	SNR	ω_cut/(2π)	δ_resolution
(mm)	(s)	(−)	(s)	(μm)	(−)	(mHz)	(mm)
7.2	0.5	900	0.1	17	587	36	3.5
11.7	1	1100	0.2	24	930	16	5.4

z	t_p	N	Δt	Δz	SNR	ω_cut/(2π)	δ_resolution
(mm)	(s)	(−)	(s)	(μm)	(−)	(mHz)	(mm)
7.2	0.5	900	0.1	17	587	36	3.5
11.7	1	1100	0.2	24	930	16	5.4

We invert these experimental temperature signals into virtual wave signals to reconstruct the location, the width, and amplitude of the internal heat sources, showing the capability of incorporating prior information in the regularization. The reconstructed depth profile of the initial temperature using the direct regularization method T-SVD⁴ is shown in Fig. 10(b). The comparison of the different reconstructed depth profiles clearly shows the entropy production caused by thermal diffusion that is equal to information loss with the result of lower axial resolution for the deeper defect, as listed in Table III. The SNR, which was derived from the singular values, is nearly doubled for the defect depth of 11.7 mm, as the duration $t_{p}$ of the heating pulse is increased from 0.5 s to 1 s. The truncation frequency $ω_{cut}$ and the spatial resolution $δ_{resolution}$ are derived by Eqs. (28) and (29), respectively.

Using additional information about the system, e.g., through regularization constraints improves the definition of the system, decreases the available states, and leads to lower entropy, which leads to an increased spatial resolution in image reconstruction. To increase the information content in the regularization process, we introduce prior knowledge in the form of positivity, which is reasonable due to the solely non-negative amplitude values of the virtual wave, as well as for the photoacoustic pressure, in the 1D regime.⁶⁹ In contrast to the bipolar reconstruction with T-SVD, non-negative values for the reconstructed virtual wave field can be enforced using iterative regularization procedures, such as the Alternating Direction Method of Multipliers (ADMM).^70,71 Additionally, for a majority of practical problems in photothermal imaging, we can assume sparsity; since in the application of defect detection, mostly isolated singular defects occur, e.g., cracks, pores, delaminations, or inclusions of other materials, or in the application of parameter estimation, only the positions of the sources, the back wall, or individual inner boundary layers should be identified. Consequently, we have only a few point or line scatterers, which leads to a sparse virtual wave signal. Sparsity is introduced in the ADMM by an appropriate formulation of the objective function using $ℓ^{1}$ -norm minimization. As shown in the reconstructed depth profiles in Fig. 10(b), iterative regularization with prior information leads to an improvement in energy localization and hence to a higher axial resolution.

B. Active thermographic computed tomography (ATCT)

Active thermographic computed tomography (ATCT) is a new hybrid reconstruction technique that utilizes the photothermal (PT) effect,^72,73 or other thermal excitation sources as inductive heating,^25,27 mechanical friction through ultrasound, microwaves, for signal generation. In the case of photothermal computed tomography (PTCT), an optical pulse is used to irradiate an object, such as biological tissue or manufactured materials and products, to generate thermal waves within the object. The thermal waves diffuse to the sample surface, where the temperature change can be measured using infrared cameras. Due to the focal plane array of the infrared camera, up to $10^{6}$ signals can be recorded simultaneously, which is equivalent to an extremely large aperture in contrast to photoacoustic transducers. The aim of PTCT is to reconstruct an image that represents a map of the initial temperature distribution within the object from the measured photothermally induced thermal signals. The initial temperature distribution is proportional to the absorbed optical energy, which can reveal useful information of the internal structure, such as material defects, or inhomogeneous material parameters.

The reconstruction process in ATCT is a two step algorithm:

Computation of the virtual wave field: The original thermographic problem is converted to a photoacoustic imaging task by a pointwise transformation, that means separately for each pixel, of the measured temperature field into a virtual wave field as shown in Sec. IV A 2.
Virtual wave backpropagation: The resulting virtual wave field exhibits wave properties, such as wavefront propagation, reflection, and refraction. Frequency domain-synthetic aperture focusing technique (F-SAFT), a well known acoustic reconstruction method,^74,75 can be applied to image the initial temperature distribution. This allows a multidimensional image reconstruction with the advantage of a higher SNR.

The common inverse heat conduction solutions for multidimensional thermographic imaging lead to large scale optimization problems with high computation costs. The virtual wave concept for ATCT splits the severely ill-posed 3D large scale problem into a variety of small-scale 1D problems, where efficient algorithms, e.g., ADMM, can compute parallel for the solution. The 1D depth-encoded signals calculated for each location could be gathered and then reconstructed to a 3D image using F-SAFT with fast computational architecture, which could be also accelerated by graphical processing unit (GPU) programming.⁷⁶

In 3D imaging, we prefer the regularization technique ADMM including the prior information sparsity and positivity because of its improvements in terms of sensitivity and depth resolving capability compared to direct regulizers. The energy input during thermal excitation always results in positive temperature signals. By converting the temperature field into a virtual wave field, positivity cannot be applied directly for 2D or 3D wave propagation. An initially non-negative acoustic signal will take negative values during propagation.⁶⁹ The circular and spherical projections, which is in 2D the Abel transformation and in 3D the time integral of the virtual wave,¹⁸ retain positivity of the initial source. For one data point, the information gain by a positivity constraint would be only a factor of 2, but for a signal with $N$ data points this factor becomes $2^{N}$ ⁠, which can be large for higher $N$ ⁠.

In Secs. IV B 1 and IV B 2, we show two examples of ATCT with limited data (limited view), when the detectors cannot be placed around the complete object. In practical measurements, only a limited space around the specimen is available for photothermal detection. The thermal diffusivity and the resulting virtual speed of sound of the bulk are assumed to be constant and isotropic.

1. Reconstruction from a single detector plane

The experimental phantom is built up with two steel rods, which are embedded in epoxy resin (Fig. 11). The cylindrical axes of the steel rods are parallel with a distance of 5 mm and inclined to the object surface. The inductive heated steel rods work as volumetric heat sources with a heating time of $t_{p} = 2 s$ similar to the 1D photothermal example (Sec. IV A 2). The temperature field is measured on a single surface at $z = 0$ ⁠. The measured signal decreases with increasing depth of the rods due to the homogeneous diffusion in each direction. The spatial- and time discretizations were $Δ_{x} = 168 μ m$ and $Δ_{t} = 40 ms$ ⁠.

FIG. 11.

View large Download slide

Experimental setup: (a) Test specimen built up with steel rods with a diameter of 1.5 mm that are embedded in epoxy resin and in (b) the principle sketch of the measurement setup: the steel rods are stimulated by induction of eddy current. The resulting change of the surface temperature is measured with an infrared camera.

Reconstruction results based on this phantom and T-SVD regularization were also published by Burgholzer et al.²⁷ In this Tutorial, not only T-SVD is used for regularization, but also the iterative regularization scheme ADMM is used, which enables the incorporation of prior information in the form of positivity and sparsity for defect reconstruction. A cross section of the virtual wave field reconstructed by T-SVD and ADMM is illustrated in Fig. 12. The virtual wave fields were normalized to the maximum wave amplitude that was obtained by ADMM regularization. Each cross section shows the typical scattering hyperbola as a consequence of the cylindrical steel rods, but the ADMM regularization leads to a much sharper localization. Based on the virtual wave field, the initial temperature distribution can be calculated by applying inverse wave propagation methods like F-SAFT.

FIG. 12.

View large Download slide

2D cross section of the reconstructed virtual wave field of the two steel bars at $x = 3 mm$ produced by the use of two different regularization techniques, where (a) the direct method of T-SVD and (b) the iterative procedure of ADMM with prior information are used. The virtual scattering hyperbolas are visible. The reconstructed (rec.) temperature cross sections were normalized by the maximum temperature obtained with ADMM regularization (norm. units).

Figure 13 illustrates the resulting 3D initial temperature distribution for the different regularization methods. The amplitude of the reconstructed initial temperature field was binarized with a threshold value to represent the localization of thermal sources at a point within the 3D sample volume. The circles indicate the real position of the steel rods in the epoxy cylinder. Figure 13(a) exhibits the T-SVD based reconstruction. It is visible that the steel rods can be distinguished until $z = 4 mm$ ⁠, as expected for direct reconstructions without any prior information. The resolution is in the range of the depth of the sources. Moreover, artifacts are present in regions with low SNR. Figure 13(b) illustrates the reconstruction employing ADMM. In this process, positivity and sparsity were introduced to improve the regularized solution. Here, the steel rods are separable in the domain $10 < x < 40 mm$ ⁠. Moreover, the artifacts are significantly reduced compared to T-SVD.

FIG. 13.

View large Download slide

Isosurface illustration of the reconstructed internal heat sources obtained via (a) T-SVD and (b) ADMM employing F-SAFT.

2. 3D plane detection

An example of ATCT with temperature data from three orthogonal detection planes is shown in Fig. 14. The phantom is an epoxy cube with an edge length of $a = 15 mm$ containing four spherical steel spheres with a diameter of 2 mm. The epoxy material, the steel beads, the inductive heating device, and the position of the sample inside the coil were the same as in Sec. IV B 1. The spatial and temporal discretization of the IR camera were $Δ x = Δ y = Δ z = 172 μ m$ and $Δ t = 50 ms$ ⁠. The bottom and the sidewalls of the cube were thermally isolated to ensure adiabatic boundary conditions. The time-dependent temperature fields were measured on the right side plane (⁠ $x = a$ ⁠), front plane (⁠ $y = 0$ ⁠), and on the top plane (⁠ $z = a$ ⁠) for a time duration of 120 s after the pulse excitation. The pulse time of 7.5 s was chosen to obtain also thermographic signals from the deeper spheres.

FIG. 14.

View large Download slide

(a) Experimental setup and detection planes of the ATCT measurements and (b) the isosurface illustration of the reconstructed internal heat sources obtained with ADMM for three different detector planes and a superposition of the single detector plane reconstructions. The steel spheres had depths of 4.3 mm, 7.5 mm, and 10.7 mm from the detection planes.

Reconstructions based on this phantom and T-SVD regularization were also published in Burgholzer et al.²⁵ In this Tutorial and in contrast to prior works, we use ADMM regularization and incorporate prior information in the form of positivity and sparsity to improve the quality of the regularized solution. The three corresponding reconstructions using ADMM for regularization and F-SAFT for image reconstruction are shown in Fig. 14. For all three reconstruction planes, limited view artifacts are visible and the inclusions far from the detector planes cannot be reconstructed at all.

V. SUMMARY, CONCLUSIONS, AND OUTLOOK

In nature, real processes are always irreversible and show a clear direction of the arrow of time. Some processes, such as the propagation of an acoustic wave, especially for low frequencies in a medium with low viscosity, are in a good approximation reversible. The corresponding wave equation for pressure is symmetric in time; for every pressure wave, which is a solution of this equation, the time-reversed pressure wave is also a solution. Even if some acoustic waves can be described in a good approximation by the wave equation, acoustic attenuation cannot be avoided totally. Therefore, we call this ideal solution a “virtual” wave, in comparison to the real measured signal. In this Tutorial, it is demonstrated for thermal and acoustic waves that the measured and the virtual signal are linked by a local relation [Eq. (1)].

For non-destructive imaging of subsurface structures, from an information theoretical viewpoint, it is advantageous to use “nearly reversible” waves to transport the information from the structures to the sample surface, such as acoustic waves with low attenuation. Due to heat diffusion, thermal waves exhibit a higher entropy production, which is shown to be equal to the information loss (Sec. III) and thus results in a lower spatial resolution (Tables I and II or Figs. 2 and 3). In Sec. III, two simple examples for stochastic processes are given, where the close connection between fluctuation and dissipation (“kicked Brownian particle”) and diffusion (“ideal gas diffusion between two boxes”) can be shown explicitly.

In Sec. IV, the information theoretical cutoff criterion using Chernoff–Stein’s Lemma (Sec. III) is applied in the frequency domain to photoacoustic and photothermal imaging. The information related criterion gives the same truncation value as when the signal amplitude gets less than the noise level. This is always the case if the mean energy of the wave component and therefore also its information content is proportional to the square of the signal amplitude. In Sec. IV A for the 1D imaging even in fat with rather high acoustic damping, the resolution is always much better than for thermographic imaging (see also the simulation in glycerine in Sec. II A). Therefore, for thermographic imaging with this high entropy production from heat diffusion, it is essential to use additional information such as positivity and sparsity by implementing iterative algorithms, e.g., by ADMM to get a better resolution, also in the sub-mm range [Fig. 10(b)]. This has been used in thermographic computed tomography (Sec. IV B), either by 2D reconstruction from single detector plane measurements or on three perpendicular planes for a cubic sample containing steel spheres to be imaged.

Using the T-SVD without incorporating additional information for a single 1D reconstruction, the spatial resolution does not get better compared to a direct inversion of the heat diffusion. The advantage of the virtual wave concept in 1D is that more advanced regularization techniques incorporating a priori information such as sparsity or positivity can be utilized in the reconstruction process. Including additional prior information allows thermographic reconstruction with a significantly better resolution than without that additional information, as shown in Figs. 10(b), 12, and 13.^72,73 In 2D and 3D, after the computation of the virtual wave field, well known acoustic reconstruction methods such as F-SAFT are used (Sec. IV B). Here, the SNR is significantly enhanced as the heat flow in all directions is taken into account for reconstructions. In the past, one-dimensional axial image reconstruction was mainly carried out for thermographic reconstruction regardless of the lateral heat flow. The virtual wave concept respects the lateral heat flow, which yields an improved SNR and therefore an increased spatial resolution.²⁵ This is similar to averaging of many 1D measurements when imaging a layered structure, but with the essential advantage, that the virtual wave concept can be used for any 2D or 3D structure to be imaged.

For future work, deep learning neural networks are planned to be used as the second step instead of, e.g., F-SAFT. For the same data from Fig. 13, this gives another significant enhancement in resolution. In Fig. 13, the threshold level for the isosurface plot is rather critical to avoid additional artifacts. With a deep neural network trained only by simulated data, the reconstructions are very stable and choosing a threshold level for the isosurface plot is not critical as reconstruction artifacts hardly appear.⁷⁷

ACKNOWLEDGMENTS

The financial support by the Austrian Federal Ministry of Science, Research and Economy and the National Foundation for Research, Technology and Development is gratefully acknowledged. Furthermore, this work has been supported by the project “From measurement science to information gaining” by the federal government of Upper Austria, the project “Multimodal and in situ Characterization of Inhomogeneous Materials” (MiCi), by the federal government of Upper Austria and the European Regional Development Fund (EFRE) in the framework of the EU-program IWB2020. Financial support was also provided by the Austrian Research Funding Association (FFG) under the scope of the COMET programme within the research project “Photonic Sensing for Smarter Processes (PSSP)” (Contract No. 871974). This programme is promoted by BMK, BMDW, the federal state of Upper Austria, and the federal state of Styria, represented by SFG. Parts of this work have been supported by the Austrian Science Fund (FWF), Project Nos. P30747-N32 and P33019-N.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

REFERENCES

1.

C. E.

Shannon

, “

A mathematical theory of communication

,”

Bell Syst. Tech. J.

27

,

379

–

423

(

1948

).

https://doi.org/10.1002/j.1538-7305.1948.tb01338.x

Google Scholar

Crossref

2.

L.

Szilard

, “

Über die entropieverminderung in einem thermodynamischen system bei eingriffen intelligenter wesen

,”

Z. Phys.

53

,

840

–

856

(

1929

).

https://doi.org/10.1007/BF01341281

Google Scholar

Crossref

3.

R.

Landauer

, “

Information is physical

,”

Phys. Today

44

(

5

),

23

–

29

(

1991

).

https://doi.org/10.1063/1.881299

Google Scholar

Crossref

4.

P. C.

Hansen

,

Rank-Deficient and Discrete Ill-Posed Problems

(

Society for Industrial and Applied Mathematics

,

1998

).

Google Scholar

Crossref

5.

R. C.

Aster

,

B.

Borchers

, and

C. H.

Thurber

,

Parameter Estimation and Inverse Problems

, 3rd ed. (

Elsevier

,

Amsterdam

,

2018

).

Google Scholar

6.

O.

Scherzer

,

M.

Grasmair

,

H.

Grossauer

,

M.

Haltmeier

, and

F.

Lenzen

, Variational Methods in Imaging, Applied Mathematical Sciences Vol. 167 (Springer, New York, 2009).

7.

M.

Esposito

and

C.

van den Broeck

, “

Second law and Landauer principle far from equilibrium

,”

Europhys. Lett.

95

,

40004

(

2011

).

https://doi.org/10.1209/0295-5075/95/40004

Google Scholar

Crossref

8.

L.

Li

,

L.

Lin

, and

L. V.

Wang

, “

Multiscale photoacoustic tomography

,”

Opt. Photonics News

29

,

32

–

39

(

2018

).

https://doi.org/10.1364/OPN.29.4.000032

Google Scholar

Crossref

9.

P.

Burgholzer

and

G.

Hendorfer

, “

Limits of spatial resolution for thermography and other non-destructive imaging methods based on diffusion waves

,”

Int. J. Thermophys.

34

,

1617

–

1632

(

2013

).

https://doi.org/10.1007/s10765-013-1513-0

Google Scholar

Crossref

PubMed

10.

A.

Mandelis

,

Diffusion-Wave Fields

(

Springer

,

New York, NY

,

2001

).

Google Scholar

Crossref

11.

P.

Burgholzer

,

G. J.

Matt

,

M.

Haltmeier

, and

G.

Paltauf

, “

Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface

,”

Phys. Rev. E

75

,

046706

(

2007

).

https://doi.org/10.1103/PhysRevE.75.046706

Google Scholar

Crossref

12.

T. L.

Szabo

,

Diagnostic Ultrasound Imaging: Inside Out

(

Elsevier

,

San Diego

,

2014

).

Google Scholar

13.

P. J.

La Riviere

,

J.

Zhang

, and

M. A.

Anastasio

, “Image reconstruction in optoacoustic tomography accounting for frequency-dependent attenuation,” in IEEE Nuclear Science Symposium Conference Record, 23–29 October 2005 (IEEE, 2005), pp. 1841–1845.

14.

P. J.

La Rivière

,

J.

Zhang

, and

M. A.

Anastasio

, “

Image reconstruction in optoacoustic tomography for dispersive acoustic media

,”

Opt. Lett.

31

,

781

–

783

(

2006

).

https://doi.org/10.1364/OL.31.000781

Google Scholar

Crossref

PubMed

15.

H.

Ammari

,

E.

Bretin

,

V.

Jugnon

, and

A.

Wahab

, “Photoacoustic imaging for attenuating acoustic media,” in Mathematical Modeling in Biomedical Imaging II, Lecture Notes in Mathematics Vol. 2035, Mathematical Biosciences Subseries, edited by H. Ammari (Springer, Heidelberg, 2012), ISSN: 0075-8434.

16.

X. L.

Deán-Ben

,

D.

Razansky

, and

V.

Ntziachristos

, “

The effects of acoustic attenuation in optoacoustic signals

,”

Phys. Med. Biol.

56

,

6129

–

6148

(

2011

).

https://doi.org/10.1088/0031-9155/56/18/021

Google Scholar

Crossref

PubMed

17.

R.

Kowar

and

O.

Scherzer

, “Attenuation models in photoacoustics,” in Mathematical Modeling in Biomedical Imaging II, Lecture Notes in Mathematics Vol. 2035, edited by H. Ammari (Springer, Heidelberg, 2012), pp. 85–130, ISSN: 0075-8434.

18.

P.

Burgholzer

,

H.

Grün

,

M.

Haltmeier

,

R.

Nuster

, and

G.

Paltauf

, Proc. SPIE 6437, 643724 (2007).

https://doi.org/10.1117/12.700723

19.

P.

Burgholzer

,

H.

Roitner

,

J.

Bauer-Marschallinger

, and

G.

Paltauf

, Proc. SPIE 7564, 75640O (2010).

https://doi.org/10.1117/12.843220

20.

P.

Burgholzer

,

T.

Berer

,

H.

Gruen

,

H.

Roitner

,

J.

Bauer-Marschallinger

,

R.

Nuster

, and

G.

Paltauf

, “

Photoacoustic tomography using integrating line detectors

,”

J. Phys. Conf. Ser.

214

,

012009

(

2010

).

https://doi.org/10.1088/1742-6596/214/1/012009

Google Scholar

Crossref

21.

P.

Burgholzer

,

H.

Roitner

,

J.

Bauer-Marschallinger

,

H.

Grun

,

T.

Berer

, and

G.

Paltauf

, “Compensation of ultrasound attenuation in photoacoustic imaging,” in Modeling of Biological Interfacial Processes Using Thickness-Shear Mode Sensors, edited by M. G. Beghi (INTECH Open Access Publisher, 2011).

22.

B. E.

Treeby

,

E. Z.

Zhang

, and

B. T.

Cox

, “

Photoacoustic tomography in absorbing acoustic media using time reversal

,”

Inverse Probl.

26

,

115003

(

2010

).

https://doi.org/10.1088/0266-5611/26/11/115003

Google Scholar

Crossref

23.

B. E.

Treeby

and

B. T.

Cox

, “

Modeling power law absorption and dispersion for acoustic propagation using the fractional laplacian

,”

J. Acoust. Soc. Am.

127

,

2741

–

2748

(

2010

).

https://doi.org/10.1121/1.3377056

Google Scholar

Crossref

PubMed

24.

B. E.

Treeby

, “

Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering

,”

J. Biomed. Opt.

18

,

036008

(

2013

).

https://doi.org/10.1117/1.JBO.18.3.036008

Google Scholar

Crossref

PubMed

25.

P.

Burgholzer

,

M.

Thor

,

J.

Gruber

, and

G.

Mayr

, “

Three-dimensional thermographic imaging using a virtual wave concept

,”

J. Appl. Phys.

121

,

105102

(

2017

).

https://doi.org/10.1063/1.4978010

Google Scholar

Crossref

26.

P.

Burgholzer

, “

Thermodynamic limits of spatial resolution in active thermography

,”

Int. J. Thermophys.

36

,

2328

–

2341

(

2015

).

https://doi.org/10.1007/s10765-015-1890-7

Google Scholar

Crossref

PubMed

27.

P.

Burgholzer

,

G.

Stockner

, and

G.

Mayr

, “

Acoustic reconstruction for photothermal imaging

,”

Bioengineering

5

,

70

(

2018

).

https://doi.org/10.3390/bioengineering5030070

Google Scholar

Crossref

PubMed

28.

P.

Burgholzer

,

T.

Berer

,

J.

Gruber

, and

G.

Mayr

, “

Super-resolution thermographic imaging using blind structured illumination

,”

Appl. Phys. Lett.

111

,

031908

(

2017

).

https://doi.org/10.1063/1.4995410

Google Scholar

Crossref

29.

M.

Gershenson

, Proc. SPIE 10661, 106610R (2018).

https://doi.org/10.1117/12.2302712

30.

J.

Gershenson

and

M.

Gershenson

, Proc. SPIE 11004, 1100404 (2019).

https://doi.org/10.1117/12.2517415

31.

K. H.

Lee

,

G.

Liu

, and

H. F.

Morrison

, “

A new approach to modeling the electromagnetic response of conductive media

,”

Geophysics

54

,

1180

–

1192

,

https://doi.org/10.1190/1.1442753

(

1989

).

Google Scholar

Crossref

32.

K. H.

Lee

and

G.

Xie

, “

A new approach to imaging with low-frequency electromagnetic fields

,”

Geophysics

58

,

780

–

796

,

https://doi.org/10.1190/1.1443464

(

1993

).

Google Scholar

Crossref

33.

M.

Gershenson

, “

Simple interpretation of time-domain electromagnetic sounding using similarities between wave and diffusion propagation

,”

Geophysics

62

,

763

–

774

,

https://doi.org/10.1190/1.1444186

(

1997

).

Google Scholar

Crossref

34.

T. M.

Cover

and

J. A.

Thomas

,

Elements of Information Theory

, 2nd ed. (

Wiley and Chichester

,

Hoboken, NJ

,

2006

).

Google Scholar

35.

P.

Burgholzer

, “Information loss and entropy production during dissipative processes in a macroscopic system kicked out of the equilibrium” (2015), arXiv:1502.00214v2.

36.

V. E.

Gusev

and

A. A.

Karabutov

,

Laser Optoacoustics

(

American Institute of Physics

,

New York, NY

,

1993

).

Google Scholar

37.

K. P.

Kostli

,

D.

Frauchiger

,

J. J.

Niederhauser

,

G.

Paltauf

,

H. P.

Weber

, and

M.

Frenz

, “

Optoacoustic imaging using a three-dimensional reconstruction algorithm

,”

IEEE J. Sel. Top. Quantum Electron.

7

,

918

–

923

(

2001

).

https://doi.org/10.1109/2944.983294

Google Scholar

Crossref

38.

M.

Xu

and

L. V.

Wang

, “

Time-domain reconstruction for thermoacoustic tomography in a spherical geometry

,”

IEEE Trans. Med. Imaging

21

,

814

–

822

(

2002

).

https://doi.org/10.1109/TMI.2002.801176

Google Scholar

Crossref

PubMed

39.

H. S.

Carslaw

and

J. C.

Jaeger

, Conduction of Heat in Solids, 2nd ed. (Oxford Science Publications, Oxford, 1959), 1986 printing.

40.

J.

Ordonez-Miranda

,

R.

Yang

,

S.

Volz

, and

J. J.

Alvarado-Gil

, “

Steady state and modulated heat conduction in layered systems predicted by the analytical solution of the phonon Boltzmann transport equation

,”

J. Appl. Phys.

118

,

075103

(

2015

).

https://doi.org/10.1063/1.4928770

Google Scholar

Crossref

41.

Mathematical and Physical Papers, edited by G. G. Stokes (Cambridge University Press, Cambridge, 2009).

42.

V. A.

Shutilov

,

Fundamental Physics of Ultrasound

(

Gordon and Breach Science Publishers

,

New York

,

1988

).

Google Scholar

43.

S. K.

Patch

and

A.

Greenleaf

, Proc. SPIE 6437, 643726 (2007).

https://doi.org/10.1117/12.701161

44.

D.

Almond

and

P.

Patel

, Photothermal Science and Techniques, Physics and Its Applications Vol. 10 (Chapman & Hall, London, 1996).

45.

A.

Salazar

, “

Energy propagation of thermal waves

,”

Eur. J. Phys.

27

,

1349

–

1355

(

2006

).

https://doi.org/10.1088/0143-0807/27/6/009

Google Scholar

Crossref

46.

S.

Pressé

,

K.

Ghosh

,

J.

Lee

, and

K. A.

Dill

, “

Principles of maximum entropy and maximum caliber in statistical physics

,”

Rev. Mod. Phys.

85

,

1115

–

1141

(

2013

).

https://doi.org/10.1103/RevModPhys.85.1115

Google Scholar

Crossref

47.

C. A.

Gearhart

, “

Boltzmann’s atom: The great debate that launched a revolution in physics

,”

Phys. Today

54

(

9

),

53

(

2001

).

https://doi.org/10.1063/1.1420513

Google Scholar

Crossref

48.

J. M. R.

Parrondo

,

C.

van den Broeck

, and

R.

Kawai

, “

Entropy production and the arrow of time

,”

New J. Phys.

11

,

073008

(

2009

).

https://doi.org/10.1088/1367-2630/11/7/073008

Google Scholar

Crossref

49.

N. G.

van Kampen

, Stochastic Processes in Physics and Chemistry, 3rd ed., North-Holland Personal Library (Elsevier, Amsterdam, 2007).

50.

C. W.

Gardiner

, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 2nd ed., Springer Series in Synergetics Vol. 13 (Springer, Berlin, 1985).

51.

H.

Risken

, The Fokker-Planck Equation: Methods of Solution and Applications, 2nd ed., Springer Series in Synergetics Vol. 18, edited by H. Risken (Springer-Verlag, Berlin, 1989), ISSN: 0172-7389.

52.

G. E.

Uhlenbeck

and

L. S.

Ornstein

, “

On the theory of the brownian motion

,”

Phys. Rev.

36

,

823

–

841

(

1930

).

https://doi.org/10.1103/PhysRev.36.823

Google Scholar

Crossref

53.

R.

Klages

,

W.

Just

, and

C.

Jarzynski

, Nonequilibrium Statistical Physics of Small Systems, Reviews of Nonlinear Dynamics and Complexity (Wiley-VCH, Weinheim, 2013).

54.

D. S.

Lemons

and

P.

Langevin

, An Introduction to Stochastic Processes in Physics: Containing “On the Theory of Brownian Motion”, edited P. Langevin, A. Gythiel, and D. S. Lemons (Johns Hopkins University Press, Baltimore, 2002).

55.

H. B.

Callen

and

R. F.

Greene

, “

On a theorem of irreversible thermodynamics

,”

Phys. Rev.

86

,

702

–

710

(

1952

).

https://doi.org/10.1103/PhysRev.86.702

Google Scholar

Crossref

56.

H. B.

Callen

and

T. A.

Welton

, “

Irreversibility and generalized noise

,”

Phys. Rev.

83

,

34

–

40

(

1951

).

https://doi.org/10.1103/PhysRev.83.34

Google Scholar

Crossref

57.

R. F.

Greene

and

H. B.

Callen

, “

On a theorem of irreversible thermodynamics. II

,”

Phys. Rev.

88

,

1387

–

1391

(

1952

).

https://doi.org/10.1103/PhysRev.88.1387

Google Scholar

Crossref

58.

J. B.

Johnson

, “

Thermal agitation of electricity in conductors

,”

Phys. Rev.

32

,

97

–

109

(

1928

).

https://doi.org/10.1103/PhysRev.32.97

Google Scholar

Crossref

59.

H.

Nyquist

, “

Thermal agitation of electric charge in conductors

,”

Phys. Rev.

32

,

110

–

113

(

1928

).

https://doi.org/10.1103/PhysRev.32.110

Google Scholar

Crossref

60.

P.

Burgholzer

,

F.

Camacho-Gonzales

,

D.

Sponseiler

,

G.

Mayer

, and

G.

Hendorfer

, Proc. SPIE 7177, 717723 (2009).

https://doi.org/10.1117/12.809074

61.

P. M.

Morse

and

K. U.

Ingard

,

Theoretical Acoustics

(

Princeton University Press

,

Princeton

,

1986

).

Google Scholar

62.

L. D.

Landau

and

E. M.

Lifshits

,

Course of Theroitical Physics

, 3rd ed. (

Pergamon Press

,

Oxford, NY

,

1980

),

Vol. 5

.

Google Scholar

63.

S. R.

de Groot

and

P.

Mazur

,

Non-Equilibrium Thermodynamics

(

Dover Publications

,

New York

,

1984

).

Google Scholar

64.

P.

Burgholzer

,

J.

Bauer-Marschallinger

, and

M.

Haltmeier

, “

Breaking the resolution limit in photoacoustic imaging using non-negativity and sparsity

,”

Photoacoustics

19

.

100191

(

2020

).

https://doi.org/10.1016/j.pacs.2020.100191

Google Scholar

Crossref

PubMed

65.

P.

Burgholzer

,

J.

Bauer-Marschallinger

,

M.

Hettich

, and

M.

Haltmeier

, Proc. SPIE 11240, 112403Y (2020).

https://doi.org/10.1117/12.2543515

66.

P.

Burgholzer

,

J.

Bauer-Marschallinger

,

B.

Reitinger

, and

T.

Berer

, “

Resolution limits in photoacoustic imaging caused by acoustic attenuation

,”

J. Imaging

5

,

13

(

2019

).

https://doi.org/10.3390/jimaging5010013

Google Scholar

Crossref

PubMed

67.

W. J.

Parker

,

R. J.

Jenkins

,

C. P.

Butler

, and

G. L.

Abbott

, “

Flash method of determining thermal diffusivity, heat capacity, and thermal conductivity

,”

J. Appl. Phys.

32

,

1679

–

1684

(

1961

).

https://doi.org/10.1063/1.1728417

Google Scholar

Crossref

68.

S.

Breitwieser

,

G.

Zauner

, and

G.

Mayr

, “

Characterization of mid-wavelength quantum infrared cameras using the photon transfer technique

,”

Infrared Phys. Technol.

106

,

103283

(

2020

).

https://doi.org/10.1016/j.infrared.2020.103283

Google Scholar

Crossref

69.

L. V.

Wang

, Photoacoustic Imaging and Spectroscopy, Optical Science and Engineering Vol. 144 (CRC Press, Boca Raton, 2009).

70.

N.

Parikh

, “

Proximal algorithms

,”

Found. Trends Optim.

1

,

127

–

239

(

2014

).

https://doi.org/10.1561/2400000003

Google Scholar

Crossref

71.

S.

Boyd

, “

Distributed optimization and statistical learning via the alternating direction method of multipliers

,”

Found. Trends Mach. Learn.

3

,

1

–

122

(

2010

).

https://doi.org/10.1561/2200000016

Google Scholar

Crossref

72.

G.

Thummerer

,

G.

Mayr

,

P. D.

Hirsch

,

M.

Ziegler

, and

P.

Burgholzer

, “

Photothermal image reconstruction in opaque media with virtual wave backpropagation

,”

NDT E Int.

112

,

102239

(

2020

).

https://doi.org/10.1016/j.ndteint.2020.102239

Google Scholar

Crossref

73.

G.

Thummerer

,

G.

Mayr

,

M.

Haltmeier

, and

P.

Burgholzer

, “

Photoacoustic reconstruction from photothermal measurements including prior information

,”

Photoacoustics

19

,

100175

(

2020

).

https://doi.org/10.1016/j.pacs.2020.100175

Google Scholar

Crossref

PubMed

74.

D.

Lévesque

,

A.

Blouin

,

C.

Néron

, and

J.-P.

Monchalin

, “

Performance of laser-ultrasonic F-SAFT imaging

,”

Ultrasonics

40

,

1057

–

1063

(

2002

).

https://doi.org/10.1016/S0041-624X(02)00256-1

Google Scholar

Crossref

PubMed

75.

L. J.

Busse

, “

Three-dimensional imaging using a frequency-domain synthetic aperture focusing technique

,”

IEEE Trans. Ultrason. Ferroelectr. Freq. Control

39

,

174

–

179

(

1992

).

https://doi.org/10.1109/58.139112

Google Scholar

Crossref

PubMed

76.

S.

Liu

,

X.

Feng

,

F.

Gao

,

H.

Jin

,

R.

Zhang

,

Y.

Luo

, and

Y.

Zheng

, “

GPU-accelerated two dimensional synthetic aperture focusing for photoacoustic microscopy

,”

APL Photonics

3

,

026101

(

2018

).

https://doi.org/10.1063/1.5005145

Google Scholar

Crossref

77.

P.

Kovacs

,

B.

Lehner

,

G.

Thummerer

,

G.

Mayr

,

P.

Burgholzer

, and

M.

Huemer

, “A hybrid approach for thermographic imaging with deep learning,” in ICASSP 2020–2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 4 May 2020–8 May 2020 (IEEE, 2020), pp. 4277–4281.

2020

Author(s)

Linking information theory and thermodynamics to spatial resolution in photothermal and photoacoustic imaging

I. INTRODUCTION

II. PHOTOACOUSTIC AND PHOTOTHERMAL IMAGING MODEL

A. Compensation of acoustic attenuation for 1D photoacoustic imaging

B. Virtual wave concept for thermographic reconstruction

III. THERMODYNAMICS AND INFORMATION

A. Kicked Brownian particle

B. Ideal gas diffusion between two boxes

IV. IMAGING AS AN INVERSE PROBLEM

A. 1D imaging

1. Acoustic attenuation in a fat tissue layer

2. Axial thermal profiling of planar heat sources

B. Active thermographic computed tomography (ATCT)

1. Reconstruction from a single detector plane

2. 3D plane detection

V. SUMMARY, CONCLUSIONS, AND OUTLOOK

ACKNOWLEDGMENTS

DATA AVAILABILITY

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Linking information theory and thermodynamics to spatial resolution in photothermal and photoacoustic imaging

I. INTRODUCTION

II. PHOTOACOUSTIC AND PHOTOTHERMAL IMAGING MODEL

A. Compensation of acoustic attenuation for 1D photoacoustic imaging

B. Virtual wave concept for thermographic reconstruction

III. THERMODYNAMICS AND INFORMATION

A. Kicked Brownian particle

B. Ideal gas diffusion between two boxes

IV. IMAGING AS AN INVERSE PROBLEM

A. 1D imaging

1. Acoustic attenuation in a fat tissue layer

2. Axial thermal profiling of planar heat sources

B. Active thermographic computed tomography (ATCT)

1. Reconstruction from a single detector plane

2. 3D plane detection

V. SUMMARY, CONCLUSIONS, AND OUTLOOK

ACKNOWLEDGMENTS

DATA AVAILABILITY

REFERENCES

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only