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 |
Iterative Image Improvement Using Time Reversal Example
On this page… |
---|
Reconstruction using data from a line array |
Overview
This example demonstrates how photoacoustic image reconstruction may be improved iteratively. First, the sensor data is simulated using kspaceFirstOrder2D
. Then, in the first reconstruction, an image is formed from data recorded on a line array using time reversal. In the second, an image is formed from data recorded on an L-shaped sensor array (also using time reversal). In the final stage, this image is improved iteratively. This example builds on the 2D Time Reversal Reconstruction For A Line Sensor Example.
Simulating the sensor data
The example begins by defining an initial pressure distribution using a pre-defined image of the word 'k-Wave'.
% load an image for the initial pressure distribution p0_image = loadImage('k-Wave.png'); % make it binary p0_image = (p0_image>0); % smooth and scale the initial pressure distribution p0_magnitude = 2; p0 = p0_magnitude*smooth(kgrid, p0_image, true); % assign to the source structure source.p0 = p0;
The acoustic pressure time series that would be measured at each of the sensors in an L-shaped sensor array are then simulated using kspaceFirstOrder2D
.
Reconstruction using data from a line array
Part of this data - for those sensors forming a line array - is then used in a time reversal image reconstruction, just as in the 2D Time Reversal Reconstruction For A Line Sensor Example.
% define a line array of sensors sensor.mask = zeros(Nx, Ny); sensor.mask(1, :) = 1; % assign the time reversal data sensor.time_reversal_boundary_data = sensor_data((sensor.mask(original_sensor_indices) == 1),:); % set the initial pressure to be zero source.p0 = 0; % run the time reversal reconstruction p0_line = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
The initial pressure distribution and image resulting from this first reconstruction are shown below:
Reconstruction using data from an L-shaped array
The image reconstruction is then repeated using all the data measured over the L-shaped array.
% extend the array to the original L-shaped form sensor.mask(:, 1) = 1; % assign the time reversal data sensor.time_reversal_boundary_data = sensor_data; % run the time reversal reconstruction p0_L = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
The image is noticeably improved, but can be improved further still by iterating, which is described below.
Iteratively improving the image reconstruction
Beginning with the image obtained from the L-shaped sensor, an iterative procedure is now used to improve it. This improvement is achieved by adding an additional set of (imaginary) sensors so that the sensor array encloses the region of interest. The missing data for these 'sensors' is then estimated using the forward model with the best image obtained so far. (In fact, a non-negativity condition is applied to the image before it is used in the forward model, as we know that in photoacoustics the initial pressure distribution cannot be negative.) The original data and the estimated data are then used together to form an image, which - at least within the region satisfying the visibility (audibility) condition - will be an improvement on the previous image. The estimate of the missing data and subsequent image reconstruction may then be iterated for even greater image improvement.
Niterations = 5; for loop = 1:Niterations % apply a non-negativity condition source.p0 = p0_iterated.*(p0_iterated>=0); % set the sensor mask to be just the missing sides of the box sensor.mask = zeros(Nx, Ny); sensor.mask(end, :) = 1; sensor.mask(:,end) = 1; %remove the time-reversed field so the forward model will run sensor = rmfield(sensor,'time_reversal_boundary_data'); % estimate the data on the missing sides of the box estimated_sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % define the full box array sensor.mask = zeros(Nx, Ny); sensor.mask(1,:) = 1; sensor.mask(:,1) = 1; sensor.mask(end,:) = 1; sensor.mask(:,end) = 1; %combine the measured data with the newly estimated data combined_sensor_data = zeros(sum(sensor.mask(:)==1),kgrid.Nt); combined_sensor_data(sensor_box_indices1,:) = estimated_sensor_data; combined_sensor_data(sensor_box_indices2,:) = sensor_data; % assign the time reversal data sensor.time_reversal_boundary_data = combined_sensor_data; % reset the initial pressure source.p0 = 0; % run the time reversal reconstruction p0_iterated = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); end
The images from the non-iterated and iterated reconstructions using the data from the L-shaped array are shown above. Notice that within the region in which the visibility condition is satisfied (ie. the normals to every boundary in the image within this region cross the measurement surface somewhere) the image is recovered well, but outside this region there remain artefacts. That is inevitable and cannot be improved without recording more data (or perhaps making other a priori assumptions).
To show that the accuracy of the amplitudes is significantly improved, a horizontal profile through the centre of the image is shown below for the various reconstructions.
Image Reconstruction With Bandlimited Sensors | Attenuation Compensation Using Time Reversal |
© 2009-2014 Bradley Treeby and Ben Cox.