k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
- Getting Started
- Examples
- Initial Value Problems
- Example: Homogenous Propagation Medium
- Example: Using A Binary Sensor Mask
- Example: Defining A Sensor Mask By Opposing Corners
- Example: Loading External Image Maps
- Example: Heterogeneous Propagation Medium
- Example: Saving Movie Files
- Example: Recording The Particle Velocity
- Example: Defining A Gaussian Sensor Frequency Response
- Example: Comparison Of Modelling Functions
- Example: Setting An Initial Pressure Gradient
- Example: Simulations In One Dimension
- Example: Simulations In Three Dimensions
- Example: Photoacoustic Waveforms in 1D, 2D and 3D
- Time Varying Source Problems
- Example: Monopole Point Source In A Homogeneous Propagation Medium
- Example: Dipole Point Source In A Homogeneous Propagation Medium
- Example: Simulating Transducer Field Patterns
- Example: Steering A Linear Array
- Example: Snell's Law And Critical Angle Reflection
- Example: The Doppler Effect
- Example: Diffraction Through A Slit
- Example: Simulations In Three-Dimensions
- Sensor Directivity
- Example: Focussed Detector in 2D
- Example: Focussed Detector in 3D
- Example: Modelling Sensor Directivity in 2D
- Example: Modelling Sensor Directivity in 3D
- Example: Sensor Element Directivity in 2D
- Example: Focussed 2D Array with Directional Elements
- Photoacoustic Image Reconstruction
- Example: 2D FFT Reconstruction For A Line Sensor
- Example: 3D FFT Reconstruction For A Planar Sensor
- Example: 2D Time Reversal For A Line Sensor
- Example: 2D Time Reversal For A Circular Sensor
- Example: 3D Time Reversal For A Planar Sensor
- Example: 3D Time Reversal For A Spherical Sensor
- Example: Image Reconstruction With Directional Sensors
- Example: Image Reconstruction With Bandlimited Sensors
- Example: Iterative Image Improvement Using Time Reversal
- Example: Attenuation Compensation Using Time Reversal
- Example: Attenuation Compensation Using Time Variant Filtering
- Example: Automatic Sound Speed Selection
- Diagnostic Ultrasound Simulation
- Example: Defining An Ultrasound Transducer
- Example: Simulating Ultrasound Beam Patterns
- Example: Using An Ultrasound Transducer As A Sensor
- Example: Simulating B-mode Ultrasound Images
- Example: Simulating B-mode Images Using A Phased Array
- Numerical Analysis
- Example: Controlling The Absorbing Boundary Layer
- Example: Source Smoothing
- Example: Filtering A Delta Function Input Signal
- Example: Modelling Power Law Absorption
- Example: Modelling Nonlinear Wave Propagation
- Example: Optimising k-Wave Performance
- Using The C++ Code
- Elastic Wave Propagation
- Example: Explosive Source In A Layered Medium
- Example: Plane Wave Absorption
- Example: Shear Waves And Critical Angle Reflection
- Example: Simulations In Three Dimensions
- Functions - By Category
- Functions - Alphabetical List
- Release Notes
- License
k-Wave Toolbox |
2D FFT Reconstruction For A Line Sensor Example
On this page… |
---|
Overview
This example demonstrates the use of k-Wave for the reconstruction of a two-dimensional photoacoustic wave-field recorded over a linear array of sensor elements. The sensor data is simulated using kspaceFirstOrder2D
and reconstructed using kspaceLineRecon
. It builds on the Homogeneous Propagation Medium and Heterogeneous Propagation Medium examples.
Defining a smooth initial pressure distribution
The sensor data used for the reconstruction is generated using kspaceFirstOrder2D
in the same way as in the Initial Value Problems examples. (The reconstruction of experimental data can approached in the same fashion by substituting the numerical sensor data with experimental measurements.) It is important to note that, by default, the function kspaceFirstOrder2D
spatially smooths the initial pressure distribution using smooth
(see the Source Smoothing Example). To provide a valid comparison, the output of the reconstruction code kspaceLineRecon
must be compared to the smoothed version of the initial pressure used in the simulation. This can be achieved by explicitly using smooth
before passing the initial pressure to kspaceFirstOrder2D
.
% smooth the initial pressure distribution and restore the magnitude source.p0 = smooth(kgrid, source.p0, true);
The default smoothing within kspaceFirstOrder2D
can be modified via the optional input parameter 'Smooth'
. If given as a single Boolean value, this setting controls the smoothing of the initial pressure along with the sound speed and density distributions. To control the smoothing of these distributions individually, 'Smooth'
can also be given as a three element array controlling the smoothing of the initial pressure, the sound speed, and the density, respectively.
Defining the time array
Although kspaceFirstOrder2D
can be used to automatically define the simulation length and time-step (by default kgrid.t_array
is set to 'auto'
by makeGrid
), this information is also required by kspaceLineRecon
and thus in this example kgrid.t_array
must be explicitly created. This can be easily achieved by directly calling makeTime
. This is the same function that is called internally by the first-order simulation functions when kgrid.t_array
is set to 'auto'
.
% create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
Simulating the sensor data
The FFT reconstruction function kspaceLineRecon
requires data recorded along an equally spaced line-shaped array of sensor points. A sensor with this shape can be created by defining a binary sensor mask matrix with a line of 1's along the first matrix row (i.e., a line along x = const).
% define a binary line sensor sensor.mask = zeros(Nx, Ny); sensor.mask(1, :) = 1;
The initial pressure distribution (along with the binary sensor mask) and the recorded sensor data returned after running the simulation are shown below.
Performing the reconstruction
The reconstruction is invoked by calling kspaceLineRecon
with the recorded sensor data, as well as the properties of the acoustic medium and the sampling parameters. By default, the sensor data input must be indexed as p(time, sensor_position)
. This means the simulated sensor data returned by kspaceFirstOrder2D
which is indexed as sensor_data(sensor_position, time)
must first be rotated. Alternatively, the optional input parameter 'DataOrder'
can be set to 'yt'
(the default settings is 'ty'
). By setting the optional input parameter 'Plot'
to true
, a plot of the reconstructed initial pressure distribution is also produced.
% reconstruct the initial pressure p_xy = kspaceLineRecon(sensor_data.', dy, dt, medium.sound_speed, 'Plot', true);
As the reconstruction runs, the size of the recorded data and the time to compute the reconstruction are both printed to the command line.
Running k-Wave line reconstruction... grid size: 216 by 778 grid points interpolation mode: *nearest computation completed in 0.11074s
Regardless of the physical alignment of the sensor within the acoustic medium, the reconstruction is always returned as if the sensor was located along the first matrix row (i.e., x = const). The resolution of the reconstruction in the y-direction is defined by the physical location and spacing of the sensor elements while the resolution in the x-direction is defined by the sampling rate at which the pressure field is recorded (i.e., dt). Consequently, the reconstructed initial pressure map typically will have a much finer discretisation in the x- (time) direction. By default, the reconstructed initial pressure distribution is not re-scaled or thresholded in any way. However, a positivity condition can be automatically enforced by setting the optional input parameter 'PosCond'
to true
.
The limited view problem
To directly compare the initial pressure distribution with that produced from the reconstruction, it is convenient to interpolate the reconstructed pressure onto a k-space grid with the same dimensions as source.p0
.
% define a second k-space grid using the dimensions of p_xy [Nx_recon, Ny_recon] = size(p_xy); kgrid_recon = makeGrid(Nx_recon, dt*medium.sound_speed, Ny_recon, dy); % resample p_xy to be the same size as source.p0 p_xy_rs = interp2(kgrid_recon.y, kgrid_recon.x - min(kgrid_recon.x(:)), p_xy, kgrid.y, kgrid.x - min(kgrid.x(:)));
The interpolated reconstructed initial pressure distribution (using a positivity condition) with the same plot scaling as the plot of source.p0
above is shown below. A x = const slice through the center of the larger disc is also shown for comparison. The reconstructed pressure magnitude is decreased due to the limited view problem; the reconstruction is only exact if the data is collected over an infinite line, here a finite length sensor is used.
Setting the interpolation mode
The FFT reconstruction algorithm relies on the interpolation between a temporal and a spatial domain coordinate with different inherent spacings. Both the speed and accuracy of the reconstruction is thus dependent on the method used for this interpolation. This can be controlled via the optional input parameter 'Interp'
which is passed directly to interp2
. By default, this is set to '*nearest'
which optimises the interpolation for speed. The accuracy of the interpolation can be improved by setting 'Interp'
to '*linear'
or '*cubic'
, and re-running the simulation. This increases the time to compute the reconstruction, however, the noise in the reconstruction is noticeably improved, and the error in the magnitude of the reconstruction is also slightly reduced.
In practice, '*linear'
interpolation provides a good balance between reconstruction speed and image artefacts.
Photoacoustic Image Reconstruction | 3D FFT Reconstruction For A Planar Sensor |
© 2009-2014 Bradley Treeby and Ben Cox.