k-Wave Toolbox |
![]() ![]() |
On this page… |
---|
First order acoustic equations Initial value problem for photoacoustics |
When modelling linear sound propagation in soft biological tissue, it is usually assumed that the propagation medium is isotropic, there is no net flow and that shear waves can be neglected. Under these conditions the acoustic pressure, p, acoustic particle velocity, u, and acoustic density, rho, are related by three first order equations corresponding to momentum conservation (sometimes called the equation of motion or Euler's equation but essentially Newton's second law), mass conservation (sometimes called the continuity equation), and an equation of state:
These can be rearranged to give the linear acoustic wave equation:
where the sound speed, c, and ambient density, rho_0, can vary with position x, and the pressure, p, will depend in general on both position and time, p(x,t).
k-Wave solves the coupled first order acoustic equations using a k-space approach, described below. The equations can be solved either in conjunction with a space- and time-varying pressure source term, or with an initial condition for the pressure. The latter is a good model for photoacoustic wave generation.
The forward problem in photoacoustics is to calculate the acoustic (ultrasonic) field as a function of time following the absorption of an optical pulse by an elastic fluid with heterogeneous optical absorption. (When RF radiation is used in place of optical, the term thermoacoustic is often preferred.) If, by absorbing the optical pulse, a region of fluid is heated sufficiently quickly to generate a local pressure increase, then an acoustic wave will result. When the heating occurs on a time-scale much shorter than the time necessary for the tissue to relax (for the density to change), the forward problem reduces to an initial value problem (IVP), with initial conditions
where the pressure distribution at time t=0 is called the initial (acoustic) pressure distribution. Exact solutions to this IVP are only known in simple cases, so in general it is necessary to seek a numerical solution. Typically, this has been tackled using finite difference time-domain (FDTD) algorithms. The forward model included in this toolbox uses instead a computationally efficient k-space approach, a technique closely allied to pseudo-spectral methods. It is can be used to model wave propagation in one, two, or three dimensional space.
The most commonly-used numerical methods for solving partial differential equations are finite-difference (FD) and finite-element methods. Although excellent for many applications, for time-domain modelling of high-frequency wave applications they become cumbersome and slow. This is due to the requirements for many mesh points per wavelength and small time steps to minimize unwanted numerical dispersion. The pseudo-spectral method can help reduce the first of these problems, and the k-space approach can help to overcome the second.
In a simple FD scheme, the gradient of the field is estimated using linear interpolation between its values at the mesh modes. A better estimate of the gradient could be obtained by fitting higher-order polynomials to a greater number of nodes, and calculating the derivative of the polynomial. The more points used, the higher the degree of polynomial required, and the more accurate the estimate of the derivative. The pseudo-spectral method takes this idea further and fits a Fourier series to all of the data - it is therefore sometimes referred to as a global, rather than local, method. There are two significant advantages to using Fourier series: first, the amplitudes of the Fourier components can be calculated efficiently using the Fast Fourier transform (FFT), and second, the basis functions are sinusoidal, so only two nodes per wavelength are required, rather than 10 or so in other methods.
While the pseudospectral method improves the efficiency in the spatial domain, conventional finite-difference schemes are still necessary to calculate the gradients (rates of change) in the time domain. The approximation thus incorporated introduces instability into the numerical propagation, which can only be controlled by limiting the size of the time-step. The techniques broadly classed as k-space methods attempt to relax this limitation in order to allow larger time-steps to be used without compromising accuracy. By comparing a simple pseudo-spectral time-domain model, for acoustically-homogeneous media, to an exact solution to the wave equation, it is possible to find replacement expressions for either the time-step or the spatial derivative such that the numerical solutions are exact for arbitrary large time steps. These "k-space" adjustments also lead to much improved numerical stability in the case of acoustically-heterogeneous media.
The image reconstruction problem in photoacoustic imaging is an acoustic inverse source problem: estimate the initial pressure distribution (or photoacoustic source region) given a set of acoustic pressure time histories recorded on a measurement surface. A number of algorithms have been proposed to solve this problem. Two are included here, as they are closely related to the k-space forward models described above:
In time-reversal image reconstruction, the k-space forward model (which is symmetrical with respect to time) is used with zero initial conditions, but with the measured acoustic pressure time histories imposed as a boundary condition at the position of the detectors on the measurement surface. When the measurement surface completely encloses the initial pressure distribution, this model will rebuild the acoustic pressure field inside the measurement surface in time-reversed order. When this is continued back to time t=0, the final pressure field will be the initial pressure distribution - the required photoacoustic image. Examples are given for two or three dimensional mediums.
When the measurement surface is a plane and the medium is acoustically homogeneous, a much quicker algorithm is available that calculates the initial pressure distribution in a single step. In essence, the algorithm maps the measured data, recorded as a function of position on the plane (2D) and time, to the image (3D); it maps the time-domain information into a spatial dimension. (It relates the temporal frequency to spatial frequency in the depth direction, which are linked by a dispersion relation.) The advantage of this algorithm is that it can be implemented using FFTs, and is thus very efficient. Examples are given for two or three dimensional mediums.
![]() |
Product Overview | Reference List | ![]() |
© 2009, 2010 Bradley Treeby and Ben Cox.