Increasing the number of active elements and grid points Ny

Hello,

I’m trying to simulate a B-mode image using a phantom with an inclusion. When I use 32 active transducer elements with Ny = 364, the simulation works well. Please find below the resultant B-mode image:

https://photos.app.goo.gl/8LYRwA2CR41abxcK8

However, when I change the number of active elements to 64 or 128 and also increase Ny to 652 (for 64 active elements), I get only a dark patch. I cannot find the reason, it seems the wave is not traveled properly. I’m only changing Ny and the number of active elements and not the Nx or travel time. Please find below the resultant B-mode image:

https://photos.app.goo.gl/8LYRwA2CR41abxcK8

I appreciate it if you can help me with this problem and clarify why this happens.

Please see below the codes for both images. In the second code, as I mentioned before, only Ny and the number of active elements have been changed.

Best,
Marjan

--------------------------------------------------------------------------------------------------------------------------

Code related to the first image:

% Simulating B-mode Ultrasound Images

clearvars;
clc;

% simulation settings
DATA_CAST = 'gpuArray-single'; % set to 'single' or 'gpuArray-single' to speed up computations
RUN_SIMULATION = true; % set to false to reload previous results instead of running simulation

% =========================================================================
% 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 = 1; % [grid points]

c_background = 1530; %[m/s]
c_target = 1600; %[m/s]
c0 = [c_background c_target];

f_max = 10e6; %5e6; % Maximum frequency
c0_min = min([c_background c_target]); % Mininmum apeed of sound
c0_max = max([c_background c_target]); % Maximum apeed of sound
c0_mean = mean([c_background c_target]);

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

points_per_wavelength = 5; %4;
dx = c0_min/(points_per_wavelength*f_max);
dy = dx; % [m]
dz = dx; % [m]
Nx = round(x/dx)
% Ny = round(y/dy)
Nx = 1024 - 2 * pml_x_size; % Check for Prime Factors!
% Ny = 672 - 2 * pml_y_size; % 1350 - 2 * pml_y_size; % Check for Prime Factors!

% set total number of grid points not including the PML
% Nx = 512 - 2 * pml_x_size; % [grid points]
Ny = 384 - 2 * pml_y_size; % [grid points]
Nz = 2;

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

% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================

% define the properties of the propagation medium

rho0= 1000; % [kg/m^3]

medium.alpha_coeff = 0.5; %0.73; %0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.05;
% medium.BonA = 9.21; %6;
% medium.alpha_mode = 'no_dispersion';

% create the time array
t_end = (Nx * dx) * 2.2 / c0_min; % [s]
% t_end = (Nx * dx) * 4.2 / c0_min; % [s]
kgrid.makeTime(c0, [], t_end);

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

% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = f_max; %7.5e6; %1.5e6; % [Hz]
tone_burst_cycles = 4;

fs = 1/kgrid.dt; %sampling frequency [Hz]

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

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

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

% physical properties of the transducer (considering L14-5/38 Linear
% transducer)
transducer.number_elements = 32; % total number of transducer elements
transducer.element_width = 10; % width of each element [grid points] Pitch = 0.3 mm
transducer.element_length = 1; % length of each element [grid points] Elevation height = 4 mm
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_mean; % sound speed [m/s]
transducer.focus_distance = 10e-3; % focus distance [m]
transducer.elevation_focus_distance = 16e-3; %19e-3; % 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 = 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 = kWaveTransducer(kgrid, transducer);

% print out transducer properties
transducer.properties;

% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================

% define a large image size to move across

y_FOV = 40e-3;
% number_scan_lines = 96;
number_scan_lines = round((dy+y_FOV)/(transducer.element_width*dy));
Nx_tot = Nx;
Ny_tot = Ny + number_scan_lines * transducer.element_width;
Nz_tot = Nz;

% define a random distribution of scatterers for the medium and the
% properties
background_map_mean = 1;
% background_map_std_background = 0.0026;
background_map_std_background = 0.008; %0.04;
background_map = background_map_mean + background_map_std_background * randn([Nx_tot, Ny_tot, Nz_tot]);

sound_speed_map = c_background * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;

% background_map_std_target = 0.0036;
background_map_std_target = 0.001; %0;
background_map_target = background_map_mean + background_map_std_target * randn([Nx_tot, Ny_tot, Nz_tot]);
sound_speed_target = c_target * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map_target;
density_target = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map_target;

radius_x = 5e-3;
radius_y = 5e-3;
radius_z = 0.25e-3;
r_x = round(radius_x/dx);
r_y = round(radius_y/dy);
r_z = round(radius_z/dz);

x_pos = 15e-3; % [m]
y_pos = 20e-3; % [m]
c_x = round(x_pos/dx) ;
c_y = round(y_pos/dy)+ Ny/2 + 1;
c_z = Nz_tot/2;

target_region = makeEllipsoid([Nx_tot, Ny_tot, Nz_tot], [c_x, c_y, c_z], [r_x, r_y, r_z]);
sound_speed_map(target_region == 1) = sound_speed_target(target_region == 1);
density_map(target_region == 1) = density_target(target_region == 1);

y_FOV_mm = (Ny_tot-Ny-1)*dy*1000

% =========================================================================
% RUN THE SIMULATION
% =========================================================================

% preallocate the storage
scan_lines_pr= 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, ...
'Smooth', [true,true,true]};

% run the simulation if set to true, otherwise, load previous results from
% disk
if RUN_SIMULATION

% set medium position
medium_position = 1;

%%% predeformation medium

% loop through the scan lines
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) ' (Predeformation)']);

% 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, :);

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

% extract the scan line from the sensor data
scan_lines_pre(scan_line_index, :) = transducer.scan_line(sensor_data);
save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_2.mat scan_lines_pre

% update medium position
medium_position = medium_position + transducer.element_width;

end

% save the scan lines to disk

save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_2.mat scan_lines_pre

else

% load the scan lines from disk
load Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus.mat;

end

% =========================================================================
% PROCESS THE RESULTS
% =========================================================================

% -----------------------------
% Remove Input Signal
% -----------------------------

% Precompression
% 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_pre = bsxfun(@times, scan_line_win, scan_lines_pre);

% store a copy of the middle scan line to illustrate the effects of each
% processing step

if mod(number_scan_lines,2)==0
scan_line_example_pre(1, :) = scan_lines_pre(end/2, :);
else
scan_line_example_pre(1, :) = scan_lines_pre((end+1)/2, :);
end

% -----------------------------
% 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_max * ( (1:length(kgrid.t_array)) * kgrid.dt - t0 ) / 2; % [m]

% define absorption value and convert to correct units
tgc_alpha_db_cm = medium.alpha_coeff * (tone_burst_freq * 1e-6)^medium.alpha_power;
tgc_alpha_np_m = tgc_alpha_db_cm / 8.686 * 100;

% create time gain compensation function based on attenuation value and
% round trip distance
tgc = exp(tgc_alpha_np_m * 2 * r);

% apply the time gain compensation to each of the scan lines
scan_lines_pre = bsxfun(@times, tgc, scan_lines_pre);

% store a copy of the middle scan line to illustrate the effects of each
% processing step
if mod(number_scan_lines,2)==0
scan_line_example_pre(2, :) = scan_lines_pre(end/2, :);
else
scan_line_example_pre(2, :) = scan_lines_pre((end+1)/2, :);
end

% -----------------------------
% Frequency Filtering
% -----------------------------

% filter the scan lines using both the transmit frequency and the second
% harmonic
% scan_lines_fund_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, tone_burst_freq, 100, true);
scan_lines_fund_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, tone_burst_freq, 100);
% set(gca, 'XLim', [0, 20]); %set(gca, 'XLim', [0, 6]);
% title('Predeformation');
% scan_lines_harm_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, 2 * tone_burst_freq, 30, true);
scan_lines_harm_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, 2 * tone_burst_freq, 30);
% set(gca, 'XLim', [0, 20]); %set(gca, 'XLim', [0, 6]);
% title('Predeformation');

% store a copy of the middle scan line to illustrate the effects of each
% processing step
if mod(number_scan_lines,2)==0
scan_line_example_pre(3, :) = scan_lines_fund_pre(end/2, :);
else
scan_line_example_pre(3, :) = scan_lines_fund_pre((end+1)/2, :);
end

% -----------------------------
% Envelope Detection
% -----------------------------

% envelope detection
scan_lines_fund_pre = envelopeDetection(scan_lines_fund_pre);
scan_lines_harm_pre = envelopeDetection(scan_lines_harm_pre);

% store a copy of the middle scan line to illustrate the effects of each
% processing step
if mod(number_scan_lines,2)==0
scan_line_example_pre(4, :) = scan_lines_fund_pre(end/2, :);
else
scan_line_example_pre(4, :) = scan_lines_fund_pre((end+1)/2, :);
end

% -----------------------------
% Log Compression
% -----------------------------

% normalised log compression
compression_ratio = 30;
scan_lines_fund_pre = logCompression(scan_lines_fund_pre, compression_ratio, true);
scan_lines_harm_pre = logCompression(scan_lines_harm_pre, compression_ratio, true);

% store a copy of the middle scan line to illustrate the effects of each
% processing step
if mod(number_scan_lines,2)==0
scan_line_example_pre(5, :) = scan_lines_fund_pre(end/2, :);
else
scan_line_example_pre(5, :) = scan_lines_fund_pre((end+1)/2, :);
end

% -----------------------------
% Scan Conversion
% -----------------------------

% upsample the image using linear interpolation
scale_factor = 30;
scan_lines_fund_pre = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund_pre, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
scan_lines_harm_pre = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm_pre, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');

% =========================================================================
% VISUALISATION
% =========================================================================

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)...
.* density_map(:, 1 + Ny/2:end - Ny/2, Nz/2));
axis image;
colormap(gray);
title('Impedance Image');
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_pre);
xlabel('Time [\mus]');
title('Predeformation');
set(gca, 'XLim', [5, t_end * 1e6]);

% plot the processed b-mode ultrasound image
figure;
horz_axis = (0:length(scan_lines_fund_pre(:, 1)) - 1) * transducer.element_width * dy / scale_factor * 1e3;
imagesc(horz_axis, r * 1e3, scan_lines_fund_pre.');
axis image;
colormap(gray);
set(gca, 'YLim', [2, 30]);
title('B-mode Image (Predeformation)');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Code related to the second image:

% Simulating B-mode Ultrasound Images

clearvars;
clc;

% =========================================================================
% 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 = 1; % [grid points]

c_background = 1530; %[m/s]
c_target = 1600; %[m/s]
c0 = [c_background c_target];

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

points_per_wavelength = 5; %4;
dx = c0_min/(points_per_wavelength*f_max);
dy = dx; % [m]
dz = dx; % [m]
Nx = round(x/dx)
% Ny = round(y/dy)
Nx = 1024 - 2 * pml_x_size; % Check for Prime Factors!
Ny = 672 - 2 * pml_y_size; % 1350 - 2 * pml_y_size; % Check for Prime Factors!

% set total number of grid points not including the PML
% Nx = 512 - 2 * pml_x_size; % [grid points]
% Ny = 384 - 2 * pml_y_size; % [grid points]
Nz = 2;

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

% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================

% define the properties of the propagation medium

rho0= 1000; % [kg/m^3]

medium.alpha_coeff = 0.5; %0.73; %0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.05;
% medium.BonA = 9.21; %6;
% medium.alpha_mode = 'no_dispersion';

% create the time array
t_end = (Nx * dx) * 2.2 / c0_min; % [s]
% t_end = (Nx * dx) * 4.2 / c0_min; % [s]
kgrid.makeTime(c0, [], t_end);

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

% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = f_max; %7.5e6; %1.5e6; % [Hz]
tone_burst_cycles = 4;

fs = 1/kgrid.dt; %sampling frequency [Hz]

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

% physical properties of the transducer (considering L14-5/38 Linear
% transducer)
transducer.number_elements = 64; % total number of transducer elements
transducer.element_width = 10; % width of each element [grid points] Pitch = 0.3 mm
transducer.element_length = 1; % length of each element [grid points] Elevation height = 4 mm
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]);

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

% define the transducer elements that are currently active
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 = kWaveTransducer(kgrid, transducer);

% print out transducer properties
transducer.properties;

% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================

% define a large image size to move across

y_FOV = 40e-3;
% number_scan_lines = 96;
number_scan_lines = round((dy+y_FOV)/(transducer.element_width*dy));
Nx_tot = Nx;
Ny_tot = Ny + number_scan_lines * transducer.element_width;
Nz_tot = Nz;

sound_speed_map = c_background * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;

radius_x = 5e-3;
radius_y = 5e-3;
radius_z = 0.25e-3;
r_x = round(radius_x/dx);
r_y = round(radius_y/dy);
r_z = round(radius_z/dz);

x_pos = 15e-3; % [m]
y_pos = 20e-3; % [m]
c_x = round(x_pos/dx) ;
c_y = round(y_pos/dy)+ Ny/2 + 1;
c_z = Nz_tot/2;

y_FOV_mm = (Ny_tot-Ny-1)*dy*1000

% =========================================================================
% RUN THE SIMULATION
% =========================================================================

% preallocate the storage
scan_lines_pr= zeros(number_scan_lines, kgrid.Nt);

% run the simulation if set to true, otherwise, load previous results from
% disk
if RUN_SIMULATION

% set medium position
medium_position = 1;

% loop through the scan lines
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) ' (Predeformation)']);

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

% extract the scan line from the sensor data
scan_lines_pre(scan_line_index, :) = transducer.scan_line(sensor_data);
save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_3.mat scan_lines_pre

% update medium position
medium_position = medium_position + transducer.element_width;

end

% save the scan lines to disk
save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_3.mat scan_lines_pre