Hi! I would want to change the frequency of the US source from 1.5 to 3MHz. When I change it, if i plot the signal it seems that I get a signal that is not correct. In my example I define the grid with respect to the frequency so I don't think the problem could be related to the grid size. Do you know if maybe some other parameters have to be changed when you change the frequency as for example the source strenght or the number of cycles ?
thank you for any help
Hi marlenep,
Possible things to check are whether the source frequency is still supported by the grid (I'd recommend using at least 3 PPW at the highest frequency contained in the tone burst), that you have also scaled the time step appropriately so the simulation is still stable, and that the total simulation time is still appropriate to capture the behaviour you are trying to model (e.g., scattering). The source strength or number of cycles don't need to be changed.
Brad.
thank you Brad, I increased the PPW because and seems to be better. How can i set a proper time step and total simulation time? Is still ok if I use this command?
t_end = (Nx*dx)*2.2/c0;
[kgrid.t_array delta] = makeTime(kgrid, c0, [], t_end);
I am trying to simulate an ultrasound image from a phantom with similar properties to those of the soft tissues and with two holes in the middle fill with water. I have acquired real ultrasound images using a linear transducer probe, in particular is the ultrasonic L14-5/38.
I don't have such a good result from my simulation, can you give me some advice to improve it? thank you very much for any help
This is my code ( in this code I am still using 2 PPW)
clear all
clc
close all
% simulation settings
DATA_CAST = 'gpuArray-single';
%DATA_CAST='single';
RUN_SIMULATION =true;
points_per_wavelength =2;
freq=5e6;
dx = 1482.3/(points_per_wavelength*freq);
dy=dx;
dz=dy;
Nx = round((46e-3)/dx); % phantom dimensions
Ny_tot=round((50.2e-3)/dy);
Ny=66;
Nz=round((50.2e-3)/dz);
% set the size of the perfectly matched layer (PML)
pml_x_size = 30; % [grid points]
pml_y_size = 10; % [grid points]
pml_z_size = 30; % [grid points]
% % set total number of grid points not including the PML
sc = 1;
% % create the k-space grid
kgrid = makeGrid(Nx, dx,Ny, dy, Nz, dz);
o=round(36.5e-03/dx);
p=round(36.9e-03/dx);
q=round(12.9e-03/dx);
r=round(7.5e-03/dx);
circle =makeDisc(Nx,Ny_tot,o,q,r,false)+(makeDisc(Nx,Ny_tot,o,p,r,false)); %3.65,1.29,0.75//36.5e-03 36.9e-03 7.5e-03
rect=[];
for z=1:Nz
for x=1:Nx
for y=1:Ny_tot
if (circle(x,y)==0)
rect(x,y,z)=0;
else
rect(x,y,z)=1;
end
end
end
end
figure,imshow(rect(:,:,20));
% define the properties of the propagation medium %first medium
c0 = 1540; % [m/s]water
rho0 = 1068; % [kg/m^3]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.BonA = 6;
% create the time array
t_end = (Nx*dx)*2.2/c0;
[kgrid.t_array delta] = makeTime(kgrid, c0, [], t_end);
%define properties of the input signal
source_strength = 0.25e3;%1e6; % [Pa]
tone_burst_freq = freq/sc; % [Hz]
tone_burst_cycles = 5;
% 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./(c0*rho0)).*input_signal;
transducer.input_signal = input_signal;
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
uno=round((38e-03/128)/dx); % transducer size
len=round(10e-03/dx); %transducer heigth
% physical properties of the transducer
transducer.number_elements = 32; % total number of transducer elements
transducer.element_width = uno; % width of each element [grid points]
transducer.element_length = len; % length of each element [grid points]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
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([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
% properties used to derive the beamforming delays
transducer.sound_speed = c0; % sound speed [m/s]
transducer.focus_distance = 6e-3; % focus distance [m]
transducer.elevation_focus_distance = 19e-3;% focus distance in the elevation plane [m]
transducer.steering_angle = 0;
% transducer,element_pitch=0.3e-03;% steering angle [degrees]
%transducer.steering_angle_max = 32;
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
number_active_elements = 32;
transducer.active_elements = ones(transducer.number_elements, 1);
% append input signal used to drive the transducer
transducer.input_signal = input_signal;
% create the transducer using the defined settings
transducer = makeTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
transducer.plot;
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a large image size to move across
number_scan_lines = 96;
Nx_tot = Nx;
%Ny_tot = Ny + number_scan_lines*transducer.element_width;
Nz_tot = Nz;
mu0_map = zeros(Nx, Ny_tot, Nz);
mu1_map = zeros(Nx, Ny_tot, Nz);
sigma0_map = zeros(Nx, Ny_tot, Nz);
%background
sound_speed_map = 1482.3*ones(Nx, Ny_tot, Nz);%.*background_map_s;
density_map = 994*ones(Nx, Ny_tot, Nz);%.*background_map_s;
alpha_coeff_map= ones(Nx,Ny_tot,Nz);
%water
sound_speed_map(rect==1)= 1482.3;%.*background_map(Vol==6); %giusto??
density_map(rect==1) = 994;%.*background_map(Vol==6);
alpha_coeff_map(rect==1) =0.2;
mu0_map(rect==1) = 0;%0.2;
mu1_map(rect==1) = 0;%0.4;
sigma0_map(rect==1) = 0;
%soft tissue
sound_speed_map(rect==0)= 1540;%.*background_map(Vol==6); %giusto??
density_map(rect==0) = 1068;%.*background_map(Vol==6);
alpha_coeff_map(rect==0) =0.54;
mu0_map(rect==0) = 0.53;
mu1_map(rect==0) = 0.5;
sigma0_map(rect==0) = 0;
%scattering
%generate noise maps
T0 = randn(Nx, Ny_tot, Nz);
T1 = randn(Nx, Ny_tot, Nz);
% calc noise to add
S = T0.*sigma0_map + mu0_map;
rho_noise = zeros(Nx, Ny_tot, Nz);
rho_noise(T1 <= mu1_map) = S(T1 <= mu1_map);
density_map = density_map.*(1+rho_noise);
% assign to the medium inputs
medium.sound_speed = sound_speed_map;
medium.density = density_map;
medium.alpha_coeff=alpha_coeff_map;
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% preallocate the storage
scan_lines = zeros(number_scan_lines, kgrid.Nt);
% set the input settings
input_args = {...
'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
% run the simulation if set to true, otherwise, load previous results from
% disk
if RUN_SIMULATION
% set medium position
medium_position= round((Ny_tot-(128*transducer.element_width))/2);
%medium_position = 146; %%%%%%%%%%
for scan_line_index = 1:number_scan_lines
% update the command line status
disp('');
disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]);
% load the current section of the medium
medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :);
medium.density = density_map(:, medium_position:medium_position + Ny - 1, :);
medium.alpha_coeff = alpha_coeff_map(:, medium_position:medium_position + Ny - 1, :);
% run the simulation
% sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer,'DataCast','gpuArray-single', input_args{:});
sensor_data = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer,'DataCast','gpuArray-single', input_args{:});
% extract the scan line from the sensor data
scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
% update medium position
medium_position = medium_position + transducer.element_width;
end
% save the scan lines to disk
save example_us_bmode_scan_lines scan_lines;
else
% load the scan lines from disk
load example_us_bmode_scan_lines
end
% =========================================================================
% PROCESS THE RESULTS
% =========================================================================
% -----------------------------
% Remove Input Signal
% -----------------------------
% create a window to set the first part of each scan line to zero to remove
% interference from the input signal
scan_line_win = getWin(kgrid.Nt*2, 'Tukey', 'Param', 0.05).';
scan_line_win = [zeros(1, length(input_signal)*2), scan_line_win(1:end/2 - length(input_signal)*2)];
% apply the window to each of the scan lines
scan_lines = bsxfun(@times, scan_line_win, scan_lines);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
scan_line_example(1, :) = scan_lines(end/2, :);
% -----------------------------
% Time Gain Compensation
% -----------------------------
% create radius variable assuming that t0 corresponds to the middle of the
% input signal
t0 = length(input_signal)*kgrid.dt/2;
r = c0*( (1:length(kgrid.t_array))*kgrid.dt/2 - t0); % [m]
% create time gain compensation function based on attenuation value,
% transmit frequency, and round trip distance
tgc_alpha = 0.4; % [dB/(MHz cm)]
tgc = exp(2*tgc_alpha*tone_burst_freq/1e6*r*100);
% apply the time gain compensation to each of the scan lines
scan_lines = bsxfun(@times, tgc, scan_lines);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
scan_line_example(2, :) = scan_lines(end/2, :);
% -----------------------------
% Frequency Filtering
% -----------------------------
% filter the scan lines using both the transmit frequency and the second
% harmonic
scan_lines_fund = gaussianFilter(scan_lines, 1/kgrid.dt, tone_burst_freq, 100, true);
scan_lines_harm = gaussianFilter(scan_lines, 1/kgrid.dt, 2*tone_burst_freq, 30, true);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
scan_line_example(3, :) = scan_lines_fund(end/2, :);
% -----------------------------
% Envelope Detection
% -----------------------------
% envelope detection
scan_lines_fund = envelopeDetection(scan_lines_fund);
scan_lines_harm = envelopeDetection(scan_lines_harm);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
scan_line_example(4, :) = scan_lines_fund(end/2, :);
% -----------------------------
% Log Compression
% -----------------------------
% normalised log compression
compression_ratio = 3;
scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
scan_lines_harm = logCompression(scan_lines_harm, compression_ratio, true);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
scan_line_example(5, :) = scan_lines_fund(end/2, :);
% -----------------------------
% Scan Conversion
% -----------------------------
% upsample the image using linear interpolation
scale_factor = 2;
scan_lines_fund = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
scan_lines_harm = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the medium, truncated to the field of view
figure;
imagesc((0:number_scan_lines*transducer.element_width-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, sound_speed_map(:, 1 + Ny/2:end - Ny/2, Nz/2));
axis image;
colormap(gray);
set(gca, 'YLim', [0, 40]);
title('Scattering Phantom');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');
% plot the processing steps
figure;
stackedPlot(kgrid.t_array*1e6, {'1. Beamformed Signal', '2. Time Gain Compensation', '3. Frequency Filtering', '4. Envelope Detection', '5. Log Compression'}, scan_line_example);
xlabel('Time [\mus]');
set(gca, 'XLim', [5 t_end*1e6]);
% plot the processed b-mode ultrasound image
figure;
horz_axis = (0:length(scan_lines_fund(:, 1))-1)*transducer.element_width*dy/scale_factor*1e3;
imagesc(horz_axis, r*1e3, scan_lines_fund.');
axis image;
colormap(gray);
set(gca, 'YLim', [5, 40]);
title('B-mode Image');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');
% plot the processed harmonic ultrasound image
figure;
imagesc(horz_axis, r*1e3, scan_lines_harm.');
axis image;
colormap(gray);
set(gca, 'YLim', [5, 40]);
title('Harmonic Image');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');
Hi marlenep,
You should set the total simulation time to be long enough to capture the reflections of interest. The line t_end = (Nx*dx)*2.2/c0; will set the time long enough for a single round trip the length of the domain.
Regarding the time step, I would refer to Section 2.7 in the k-Wave manual, where this is discussed in detail. In particular, the sound speed you pass to makeTime (which sets dt) and the value of medium.sound_speed_ref (which sets the value of sound speed in the k-space operator; this defaults to max(medium.sound_speed(:))) will both affect the simulation.
>> I don't have such a good result from my simulation
Can you say a bit more about what the problem is?
Brad.
Hi Bradley,
thank you very much for your advice
I share my code and the resulting scan lines matrix with you at this link:
https://imperialcollegelondon.box.com/s/vu0obiih6kn1watkojfcc8x2sheezy1o
I would like to obtain a simulated ultrasound image similar to real ultrasound image ("real_us.bmp"). I already tried playing with different parameters and the best result I obtained is the one you can find in the folder. Do you think there is a way to improve it or the limitation depends only on the mri segmentation which of course approximates the real brain tissues differentiation?
thank you very much for your time
Hello :)
I recently started studies related to reconstruction of ultrasound images and I am using kwave software (the example of the b-mode using linear transducer) to implement functions as spatial compound and pulse inversion harmonic image techniques. However, when I did the harmonic image simulations for a frequency of 3.5MHz, I got results that i think are wrong. Reading the comments in the forum and the user manual (which I'm still having difficulties to understand, because the content present there goes beyond my knowledge, and I'm still learning), I realized that there was a need to change the value of dx (mainly), but changing only dx the phantom was distorted (there were no more spheres but ellipses), even so, I could not increase the value of dy as much as I increased dx because otherwise the simulations would take too long, even using the GPU. I would like to confirm that what I understand so far is correct about the grid changes and if this distortion of the phantom really happens for different values of dx and dy?
Here's the part of my code for grid changes:
'% set total number of grid points not including the PML
Nx = 512 - 2 * pml_x_size; % [grid points]
Ny = 128 - 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 = 40e-3; % [m]
dx_old = x/(256 - 2 * pml_x_size); % [grid points]
% calculate the spacing between the grid points
dx = x / Nx; % [m]
dy = dx_old; % [m]
dz = dx_old; % [m]
'
and about the part of the code to make the scattering region:
% define a sphere for a highly scattering region
radius = 6e-3; % [m]
x_pos = 27.5e-3; % [m]
y_pos = 20.5e-3; % [m]
scattering_region1 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx_old), Nz_tot/2, round(radius/dx_old));
% assign region
sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1);
density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1);
% define a sphere for a highly scattering region
radius = 5e-3; % [m]
x_pos = 30.5e-3; % [m]
y_pos = 37e-3; % [m]
scattering_region2 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx_old), Nz_tot/2, round(radius/dx_old));
% assign region
sound_speed_map(scattering_region2 == 1) = scattering_c0(scattering_region2 == 1);
density_map(scattering_region2 == 1) = scattering_rho0(scattering_region2 == 1);
% define a sphere for a highly scattering region
radius = 4.5e-3; % [m]
x_pos = 15.5e-3; % [m]
y_pos = 30.5e-3; % [m]
scattering_region3 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx_old), Nz_tot/2, round(radius/dx_old));
'
thank you!
Amanda
Hi Amanda,
Unfortunately, the geometry functions in k-Wave assume the grid is uniform. You could write your own function to get around this without too much effort by using the Cartesian coordinates of the grid points, rather than the grid coordinates.
However, I would suggest keeping dx = dy = dz, otherwise the grid will support a different maximum frequency in each direction.
Hope that helps,
Brad.
You must log in to post.