Image Artefacts - Ultrasound Imaging

k-Wave User Forum » General Acoustics » Ultrasound Simulations

Image Artefacts - Ultrasound Imaging

(8 posts) (2 voices)

Started 9 years ago by Mithun
Latest reply from Mithun

Mithun

Member
Joined: Jun '14
Posts: 5

Dear All,

I am trying to do a simple US image simulation using K-wave [Transducer as source and sensor]. Quite new to it. I created a medium and then a ball in the middle with different speed of sound and density. I expected only the reflection from the ball. Instead I am seeing a very high signal in the beginning of time trace (May be input signal?, How to get rid?) and different reverberations at almost all depths. What may be the reason for this? I put my transducer below PML. I am pasting my script below. Can anyone please throw your views on it? It would be of great help for me to move forward.

clear all;
clc
% simulation settings
DATA_CAST = 'single';

% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================

% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 20; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]

% set total number of grid points not including the PML
Nx = 512 - 2*PML_X_SIZE; % [grid points]
Ny = 512 - 2*PML_Y_SIZE; % [grid points]
Nz = 128 - 2*PML_Z_SIZE; % [grid points]

% set desired grid size in the x-direction not including the PML
x = 30e-3; % [m]

% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = dx; % [m]
dz = dx; % [m]

% create the k-space grid
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

% Total time, 2 times for US
t_end = (x/1500)*2; % [s]
[kgrid.t_array,dt] = makeTime(kgrid, 1500, [], t_end);

% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================

% define properties of the input signal
source_strength = 5e6; % [Pa]
tone_burst_freq = 1.5e6; % [Hz]
tone_burst_cycles = 2;

% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);

% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength./(1500*1000)).*input_signal;

medium.sound_speed =1500; % [m/s]
medium.density = 1000; % [kg/m^3]
% medium.alpha_coeff = 0.0022; % [dB/(MHz^y cm)]
% medium.alpha_power = 2;
% medium.BonA = 6;

% define properties

density_map = 1000*ones(Nx_tot, Ny_tot, Nz_tot);

medium.density = density_map (:,:,:);

ball_radius = round(.1e-3/dx); % [mm]

p1 = 962*makeBall(Nx, Ny, Nz,PML_X_SIZE + round(16e-3/dx),PML_Y_SIZE + round(15e-3/dx), PML_Z_SIZE + Nz/2, ball_radius);
medium.density = medium.density + p1;

% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================

% physical properties of the transducer
transducer.number_elements = 64; % total number of transducer elements
transducer.element_width = 2; % width of each element [grid points/voxels]
transducer.element_length = 24; % length of each element [grid points/voxels]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points/voxels]
transducer.radius = inf; % radius of curvature of the transducer [m]

% calculate the width of the transducer in grid points
transducer_width = transducer.number_elements*transducer.element_width ...
+ (transducer.number_elements - 1)*transducer.element_spacing;

% use this to position the transducer in the middle of the computational grid
transducer.position = round([PML_X_SIZE+1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);

% properties used to derive the beamforming delays
transducer.sound_speed = 1500; % sound speed [m/s]
transducer.focus_distance = 15e-3; % focus distance [m]
transducer.elevation_focus_distance = 16e-3 + PML_X_SIZE; % focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]

% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';

% define the transducer elements that are currently active
transducer.active_elements = zeros(transducer.number_elements, 1);
transducer.active_elements(1:transducer.number_elements) = 1;

transducer.input_signal = input_signal;
% create the transducer using the defined settings
transducer = makeTransducer(kgrid, transducer);

% print out transducer properties
transducer.properties;

% =========================================================================
% RUN THE SIMULATION
% =========================================================================
sensor.mask=transducer.all_elements_mask;
% set the input settings
input_args = {...
'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST};

% run the simulation
sensor_data = kspaceFirstOrder3DC(kgrid, medium, transducer, transducer, input_args{:});

% extract a single scan line from the sensor data using the current
% beamforming settings
scan_line = transducer.scan_line(sensor_data);
% reshape sensor data to y, z, t

sensor_data_rs = reshape(sensor_data, transducer.number_elements*transducer.element_width, transducer.element_length, length(kgrid.t_array));

Thank You
Kind Regards,
Mithun

Posted 9 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

Hi Mithun,

There is currently a difference in the way the C++ and MATLAB codes return the data when using an object of the kWaveTransducer class as a sensor. The MATLAB code returns one time series per active transducer element, while the C++ code returns one time series per grid point that forms part of the active transducer.

If using the C++ code, before using the scan_line method, you should call:

sensor_data = transducer.combine_sensor_data(sensor_data);

This will combine the sensor data to give a single time series per active transducer element, analogous to the output from the MATLAB code.

Note that transmit apodization is not currently implemented in the C++ code, i.e., the apodization is always the same as:

transducer.transmit_apodization = 'Rectangular';

(see the discussion here).

Let me know if that doesn't solve the problem.

Brad.

Posted 9 years ago #
Mithun

Member
Joined: Jun '14
Posts: 5

Dear Brad,
Thanks for the response. Great to know about apodization in c++ code. I am not able to use transducer.combine_sensor_data(sensor_data). Its throwing an error "No appropriate method, property, or field combine_sensor_data for class kWaveTransducer". Do I have to write the method explicitly? Please help me on this.

I am not going to use scan_line method as I am interested in data from each elements separately.

Thank You,
Regards,
Mithun

Posted 9 years ago #
Mithun

Member
Joined: Jun '14
Posts: 5

Dear Brad,

I fixed it by adding the code you provided in another post to KwaveTransducer class. Thats amazing, it saved lot of time. Thank you.

My problem is that I am getting a huge initial signal in all elements. This is a phase inverted bipolar signal. Is it the input signal itself? How can I get rid of it? My transducer is below PML I guess. Even after this initial signal, I have another weak reflection also from unexpected depth. Any ideas on it looking at the code above?

Thank You,
Kind Regards,
Mithun

Posted 9 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

Hi Mithun,

This function was added as part of k-Wave V1.1 (released a couple of weeks ago).

The signal at the beginning will be your input signal. Inside the code, the source and sensor act independently, so if they are in the same position, the sensor will also record the input signal. There's a couple of ways you could remove this:

(1) You could zero the first part of the signal, e.g., following the "Simulating B-mode Ultrasound Images Example".

(2) You could run a separate simulation in a homogeneous medium, and subtract the recorded signal from the signal you record when you have scatterers.

Are you sure it's a reflection, and not an edge wave from the edge of the transducer?

Brad.

Posted 9 years ago #
Mithun

Member
Joined: Jun '14
Posts: 5

Hi Brad,

Thank you for the valuable comments. Yes of-course, I can zero the first part of the signal. But I would like to understand whats happening actually. The high signal is coming not in the very beginning, but at around 100th sample for the middle element. It is a bipolar signal with negative cycle first and positive next. Then I thought it may be a reflection with phase shift of input signal with 1 cycle. How can I recognize an edge wave?

Combining sensor data is really useful. Thank you.

Regards,
Mithun

Posted 9 years ago #

Developer
Joined: Mar '10
Posts: 1,640

Hi Mithun,

If you plot the simulation, you'll see the wave from the edge of the aperture. It looks like a spherical wave (or cylindrical wave depending on your aperture shape), but varies in polarity on either side of planar wave front. I've created a very simple example to illustrate:

clear all;

% create the computational grid
Nx = 64;            % number of grid points in the x direction
Ny = 32;            % number of grid points in the y direction
Nz = 96;            % number of grid points in the z direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
dz = 0.1e-3;        % grid point spacing in the z direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

% define the properties of the propagation medium
medium.sound_speed = 1500;

% create time array
kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.5, 4e-6);

% create initial pressure distribution
source.p0 = zeros(Nx, Ny, Nz);
source.p0(20, 14:19, 38:59) = 1;  

% define sensor mask
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(Nx/2, Ny/2, Nz/4) = 1;
sensor.mask(3*Nx/4, Ny/2, Nz/2) = 1;

% input arguments
input_args = {'PlotPML', false, 'DataCast', 'single'};

% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

% plot the simulated sensor data
figure;
plot(sensor_data.')
ylabel('Pressure');
xlabel('Time Step');
legend('Off Axis', 'On Axis');

Hope that helps,

Brad.

Posted 9 years ago #

Mithun

Member
Joined: Jun '14
Posts: 5

Hi Brad,

That was of great help. Thank you.

Kind Regards,
Mithun

Posted 9 years ago #

RSS feed for this topic

Reply

You must log in to post.

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields

Image Artefacts - Ultrasound Imaging

Reply

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields