k-Wave User Forum » Topic: Artifacts in Plane Wave Simulations using kWaveArray

k-Wave User Forum » Topic: Artifacts in Plane Wave Simulations using kWaveArray http://www.k-wave.org/forum/topic/artifacts-in-plane-wave-simulations-using-kwavearray Support for the k-Wave MATLAB toolbox en-US Sun, 31 May 2026 15:19:23 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Artifacts in Plane Wave Simulations using kWaveArray" http://www.k-wave.org/forum/topic/artifacts-in-plane-wave-simulations-using-kwavearray#post-9011 Wed, 24 Jan 2024 11:17:15 +0000 Bradley Treeby 9011@http://www.k-wave.org/forum/ I haven't run your simulation, but the effects you're describing is what I would expect to happen given the way sources are modelled in k-Wave. Depending on the transducer design, you also see this in experimental measurements. One simple way to remove the transmitted wave from the receive data would be to run two simulations, one in a homogeneous medium, and then one with the medium properties you are trying to image, and subtract the RF data of the first from the second. nikunj101 on "Artifacts in Plane Wave Simulations using kWaveArray" http://www.k-wave.org/forum/topic/artifacts-in-plane-wave-simulations-using-kwavearray#post-8987 Wed, 13 Dec 2023 19:10:16 +0000 nikunj101 8987@http://www.k-wave.org/forum/ Hi, I am currently working on simulating RFData generated by plane wave emissions from a linear array transducer. In this process, I am encountering several artifacts in my simulations, particularly when recombining sensor data using kWaveArray. The primary issue involves the inadequate cancellation of signals from adjacent elements, especially noticeable when the transmit angle deviates from zero. This results in unwanted signals appearing in the simulated RFData. Additionally, I observe an edge effect. The edges of my simulated array produce spherical wave patterns that propagate across the array, creating an X-like artifact in the RFData. While similar phenomena occur in practical scenarios with linear arrays, the inherent directionality of actual transducer elements usually reduces the intensity of this artifact. In my simulation using kWave, even though the spherical sources are non-directional, I expect the line elements implemented via kWaveArray to exhibit some directionality, so I'm surprised I'm seeing as strong an artifact as I am. I have tried applying a roll-off apodization along the array edges to mitigate this issue, but I am seeking alternative solutions. Could you suggest other methods to reduce the edge effect artifact? Additionally, how might I address the interference artifact arising from neighboring elements? Here is a minimal working example of my code using a homogenous medium. You can adjust the transmit angle by changing the TXangle parameter in the first section: <code></code>` % Coherent Plane Wave Compounding Simulation Using K-Wave % % This script demonstrates how to setup a linear array transducer using the % kWaveArray class for coherent plane wave compounding in ultrasound imaging. % It involves creating a simulation grid, defining a medium with point targets, % initializing a transducer array, generating a focused source, simulating wave % propagation, and collecting sensor data. % % Author: Nikunj Khetan % Last Update: 21st November 2023 clearvars; % close all %% DEFINE LITERALS - Setting up parameters for the simulation % Selection of K-Wave code execution model model = 1; % Options: 1 - MATLAB CPU, 2 - MATLAB GPU, 3 - C++ code, 4 - CUDA code USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns % Medium parameters c0 = 1540; % Sound speed in the medium [m/s] rho0 = 1020; % Density of the medium [kg/m^3] % Source parameters source_f0 = (250/48)*1e6; % Frequency of the ultrasound source [Hz] source_amp = 1e6; % Amplitude of the ultrasound source [Pa] source_cycles = 3; % Number of cycles in the tone burst signal numEl = 128; % Number of elements in the transducer array element_width = 2.3e-4; % Width of each transducer element [m] element_pitch = 2.3e-4; % Pitch - distance between the centers of adjacent elements [m] RF_fs = source_f0*4; % Sampling Frequency of final RFData % Define transmission angles for plane wave compounding na = 1; % Number of angles for transmission TXangle = 10*pi/180; % Grid parameters grid_size_x = 40e-3; % Grid size in x-direction [m] grid_size_y = 10e-3; % Grid size in y-direction [m] % Computational parameters ppw = 4; % Points per wavelength t_end = round(grid_size_y/c0,6); % Simulation duration [s] cfl = 0.5; % Courant-Friedrichs-Lewy (CFL) number for stability %% GRID - Creating the computational grid % Calculate grid spacing based on PPW and source frequency dx = c0 / (ppw * source_f0); % Grid spacing [m] dy = dx; % Compute grid size Nx = roundEven(grid_size_x / dx); % Number of grid points in x-direction Ny = roundEven(grid_size_y / dy); % Number of grid points in y-direction % Create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dy); % Create the time array kgrid.makeTime(c0, cfl, t_end); dsFactor = (1/kgrid.dt)/RF_fs; %% SOURCE - Initializing the transducer array and source signal [karray, ElemPos] = initArray(kgrid, numEl, element_pitch, element_width); % Create source signal using a tone burst source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles); %% MEDIUM - Defining the medium properties medium.sound_speed = c0 * ones([Nx, Ny]); % sound speed [m/s] medium.density = rho0 * ones([Nx, Ny]); % density [kg/m3] %% SENSOR - Setting up the sensor mask and properties % Assign binary mask from karray to the sensor mask sensor.mask = karray.getArrayBinaryMask(kgrid); % set the record mode such that only the rms and peak values are stored sensor.record = {'p'}; % Define frequency response of the sensor sensor.frequency_response = [source_f0, 100]; %% SIMULATION - Running the simulation for different transmission angles % Preallocate arrays for time delays and RF data time_delays = zeros(na, numEl); % Simulation input options input_args = {'PMLSize', 'auto', 'PMLInside', false, 'PlotPML', false, 'DisplayMask', 'off'}; RFData = zeros(numEl, kgrid.Nt, na); % Loop over each angle for plane wave compounding for i = 1:na [source, time_delays(i, :)] = genSource(kgrid, source_f0, source_cycles, source_amp, TXangle(i), karray, ElemPos, c0); sensor_data = runSim(kgrid, medium, source, sensor, input_args, model,source_amp); RFData(:, :, i) = karray.combineSensorData(kgrid, sensor_data.p); end % Rearrange RF data dimensions for further processing RFData = downsample(flip(permute(RFData, [2, 1, 3]),2),dsFactor); figure colormap gray imagesc(log10(abs(RFData(:,:,1)))) %% HELPER FUNCTIONS function [karray, ElemPos] = initArray(kgrid, element_num, element_pitch, element_width) % Initializes the transducer array. % Args: % kgrid: The k-Wave grid object. % element_num: Number of elements in the array. % element_pitch: Distance between the centers of adjacent elements. % element_width: Width of each element. % Returns: % karray: The k-Wave array object. % ElemPos: The positions of the elements in the array. % Create empty kWaveArray object with specified BLI tolerance and upsampling rate karray = kWaveArray('BLITolerance', 0.05, 'UpsamplingRate', 10); % Calculate the center position for the first element L = element_num * element_pitch / 2; ElemPos = -(L - element_pitch / 2) + (0:element_num - 1) * element_pitch; % Add rectangular elements to the array for ind = 1:element_num % Set element position x_pos = ElemPos(ind); % Define start and end points of the element start_point = [x_pos - element_width / 2, kgrid.y_vec(1)]; end_point = [x_pos + element_width / 2, kgrid.y_vec(1)]; % Add line element to the array karray.addLineElement(start_point, end_point); end end function [source, time_delays] = genSource(kgrid, source_f0, source_cycles, source_amp, theta, karray, ElemPos, c0) % Generates the source signal with time delays for each transducer element. % Args: % kgrid: The k-Wave grid object. % source_f0: Frequency of the source. % source_cycles: Number of cycles in the tone burst signal. % source_amp: Amplitude of the source. % theta: Steering angle of the plane wave. % karray: The k-Wave array object. % ElemPos: The positions of the elements in the array. % c0: Speed of sound. % Returns: % source: The source object containing the mask and signals. % time_delays: Time delays applied to each element for focusing. % Calculate time delays for each element based on steering angle time_delays = ElemPos * sin(theta) / c0; time_delays = time_delays - min(time_delays); % Create time-varying source signals for each physical element source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles, 'SignalOffset', round(time_delays / kgrid.dt)); % Assign binary mask for the source source.p_mask = karray.getArrayBinaryMask(kgrid); % Assign source signals to the source object source.p = karray.getDistributedSourceSignal(kgrid, source_sig); end function [sensor_data] = runSim(kgrid, medium, source, sensor, input_args, model, source_amp) % Runs the simulation based on the selected model (CPU or GPU). % Args: % kgrid: The k-Wave grid object. % medium: The medium in which waves propagate. % source: The source object containing the ultrasound signal. % sensor: The sensor object to record the pressure. % input_args: Additional input arguments for the simulation. % model: The selected model for running the simulation. % Returns: % sensor_data: The recorded sensor data from the simulation. % Run the simulation based on the chosen model switch model case 1 % MATLAB CPU sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... input_args{:}, ... 'DataCast', 'single', ... 'PlotScale', [-1, 1] * source_amp); case 2 % MATLAB GPU sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... input_args{:}, ... 'DataCast', 'gpuArray-single', ... 'PlotScale', [-1, 1] * source_amp); case 3 % C++ code sensor_data = kspaceFirstOrder2DC(kgrid, medium, source, sensor, input_args{:}); case 4 % C++/CUDA GPU sensor_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}); end end <code></code>`