US probe frequency « k-Wave User Forum

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]');

k-Wave

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

US probe frequency

Reply

k-Wave

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

US probe frequency

Reply

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