I always get many NaN values in sensor_data in the linear array B mode simulation when I put two scatterers closed to each other as shown in the below figure. My simulation medium is very simple: a homogeneous medium with sound of speed as 1540 m/s and density as 1000 kg/m3 with two scatterers simulating air bubbles with sound of speed as 340 m/s and density as 1 kg/m3. When there is only one scatterer, everything is fine. However, many NaN values in the sensor data appearred when another scatterer was added to the adjacent position to the initial scatterer. If the distance between two scatterers is not zero, everything is fine, too...
Could anyone please help to solve the problem?
Below is my code:
clear
% 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 = 10; % [grid points]
pml_y_size = 5; % [grid points]
pml_z_size = 5; % [grid points]
% set total number of grid points not including the PML
Nx = 128 - 2 * pml_x_size; % [grid points]
Ny = 64 - 2 * pml_y_size; % [grid points]
Nz = 64 - 2 * pml_z_size; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 20e-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 = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
% define the properties of the propagation medium
c0 = 1540; % [m/s]
rho0 = 1000; % [kg/m^3]
medium.alpha_coeff = 0; % [dB/(MHz^y cm)]
medium.alpha_power = 0;
medium.BonA = 0;
% % create the time array
t_end = (Nx*dx)*2.5/c0; % [s]
kgrid.makeTime(c0, [], t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 1.5e6; % [Hz]
tone_burst_cycles = 4;
% 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;
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 16; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points]
transducer.element_length = 1; % 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 = 10e-3; % focus distance [m]
transducer.elevation_focus_distance = 10e-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
number_active_elements = 16;
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
number_scan_lines = 16;
Nx_tot = Nx;
Ny_tot = Ny + number_scan_lines * transducer.element_width;
Nz_tot = Nz;
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0;
background_map = background_map_mean + background_map_std * randn([Nx_tot, Ny_tot, Nz_tot]);
% define properties
sound_speed_map = c0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
simstart = 0.01;
phantom_positions = [0 0 simstart];
N_scatter = size(phantom_positions,1);
i = 1;
current_pos = phantom_positions(i,3);
x_pos = round(current_pos/dx);
sound_speed_map(x_pos, round(Ny_tot/2), round(Nz_tot/2)) = sound_speed_map(x_pos, round(Ny_tot/2), round(Nz_tot/2))*340/1540;
sound_speed_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2)) = sound_speed_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2))*340/1540;
density_map(x_pos, round(Ny_tot/2), round(Nz_tot/2)) = density_map(x_pos, round(Ny_tot/2), round(Nz_tot/2)) *0.001;
density_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2)) = density_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2)) *0.001;
% =========================================================================
% 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 = 1;
% loop through the scan lines
% for scan_line_index = 1:number_scan_lines
for scan_line_index = 1
% 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, :);
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, 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
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 - 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 = bsxfun(@times, tgc, scan_lines);
% -----------------------------
% Envelope Detection
% -----------------------------
% envelope detection
scan_lines_fund = envelopeDetection(scan_lines);
% -----------------------------
% Log Compression
% -----------------------------
% normalised log compression
compression_ratio = 3;
scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
% -----------------------------
% Scan Conversion
% -----------------------------
% upsample the image using linear interpolation
scale_factor = 1;
% =========================================================================
% 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);
colorbar;
set(gca, 'YLim', [0, 20]);
title('Scattering Phantom');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');
% 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', [0, 20]);
title('B-mode Image');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');