k-Wave Toolbox Previous   Next

Product Overview

Introduction to the k-Wave Toolbox

The k-Wave toolbox has been designed to make simulations of photoacoustic imaging simple and fast. It consists of a set of functions for solving the acoustic (ultrasonic) forward and inverse problems in photoacoustic imaging. The toolbox includes:

This beta release of the k-Wave Toolbox is focussed on modelling and image reconstruction in photoacoustic imaging, as many of the functions have been tried and tested as research tools in this area. Future releases will incorporate additional functionality of interest to the photoacoustic and wider ultrasound communities.

 Back to Top

Photoacoustic forward modelling

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. In a stationary fluid, under conditions whereby the viscosity and thermal conductivity are negligible, the acoustic pressure, p(x,t), can be modelled with the linear acoustic wave equation:

where the sound speed and density can vary with position x. The pressure p, will depend, in general, on both position x and time, t. 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.

 Back to Top

Pseudo-spectral and k-space models

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.

 Back to Top

Photoacoustic image reconstruction

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.

 Back to Top


© 2009 Bradley Treeby and Ben Cox.